Evaluation of sub-daily precipitation variability and its diurnal characteristics in satellite and reanalysis datasets against ground-based observations

The diurnal cycle of precipitation serves as a fundamental diagnostic of the climate system's energy budget and hydrological cycle. Yet comprehensive evaluations of both its regular diurnal and irregular intermittent components across multiple satellite and reanalysis products against extensive ground-based observations remain limited. This study evaluates the boreal summer diurnal cycle across five satellite products and four reanalysis datasets against 6824 ground-based stations over the Northern Hemisphere land regions, complemented by inter-product comparison over quasi-global domain between 60° S and 60° N. Analyses are performed at 3-hourly resolution, combining sub-daily variance decomposition and harmonic analysis with atmospheric moisture budget diagnosis across contrasting precipitation regimes.

Ground-based observations indicate sub-daily fluctuations account for 65.6 % of total precipitation variance, with the irregular component dominating over the repeatable diurnal signal, a fraction that all evaluated datasets systematically underestimate. Reanalysis datasets place the land precipitation maximum 3–6 h ahead of observations, driven by premature convective triggering. Among reanalysis datasets, the Japanese Reanalysis for Three Quarters of a Century (JRA-3Q) shows the smallest phase difference owing to its dynamic convective available potential energy (CAPE) generation rate (DCAPE)-based convective trigger, while the North American Regional Reanalysis (NARR) achieves superior performance over North America through direct assimilation of precipitation observations. Satellite products exhibit the observed late-afternoon precipitation maximum with a phase lag of 1–2 h, arising from the indirect detection of upper-tropospheric ice-scattering signals that persist as convective anvils expand beyond the timing of peak surface precipitation. Among satellite products, the Climate Prediction Center Morphing Technique (CMORPH) shows the closest agreement with observations, as its morphing-based propagation of passive microwave retrievals preserves precipitation timing while its probability density function (PDF)-based bias correction maintains sub-daily temporal structure. These findings provide practical guidance for dataset selection and offer actionable physical insights for improving convective parameterization schemes and satellite retrieval algorithms.

1 Introduction

Precipitation serves as the primary atmospheric engine of the global water cycle, acting as the critical link that governs the continuous exchange of moisture and latent heat between the Earth's surface and the atmosphere (Trenberth et al., 2003). While the latent heat released during condensation is a fundamental driver of global atmospheric circulation (Allan et al., 2020), precipitation itself constitutes the vital water input for terrestrial systems. As the principal source of soil moisture and groundwater, it fundamentally dictates the land-surface water budget, thereby regulating ecosystem functions and sustainable water availability (Oki and Kanae, 2006). The dynamic nature of the hydrological cycle is further underscored by precipitation extremes; excessive rainfall triggers devastating floods, while prolonged deficits lead to droughts. These deficits deplete soil moisture and shift surface energy partitioning, often culminating in complex drought–heatwave compound extremes (Miralles et al., 2019; Zscheischler et al., 2020).

The variability of precipitation spans a wide spectrum of time scales, from interannual to sub-daily, with distinct physical processes governing the hydrological cycle at each scale. Variations at interannual, seasonal, intraseasonal, and synoptic scales are associated with El Niño–Southern Oscillation (He and Li, 2019; Huang and Xie, 2015), monsoons (Wang et al., 2012; Wang and Ding, 2008), the Madden–Julian Oscillation (Jiang et al., 2020; Zhang, 2005), and extratropical cyclones (Catto et al., 2012; Hawcroft et al., 2012; Pfahl and Wernli, 2012), respectively. Beyond these large-scale oscillations, sub-daily variability constitutes a large fraction of the total precipitation budget. Satellite retrievals indicate that sub-daily internal variability explains more than 50 % of the total internal variability during the boreal summer (Ahn et al., 2022). The diurnal cycle, as the dominant component of this sub-daily variance, represents one of the most fundamental modes of the climate system, governed by the 24 h periodicity of Earth's rotation (Christopoulos and Schneider, 2021; Dai, 2001; Yang and Slingo, 2001).

Over land, the diurnal cycle of precipitation often peaks in the afternoon and is intrinsically linked to the surface energy balance. As incoming solar radiation maximizes around noon, the surface responds through intensified sensible and latent heat fluxes. The sensible heat flux deepens the planetary boundary layer (PBL), while the latent heat flux supplies the moisture to the lower atmosphere (Dai, 2024). Specifically, turbulent mixing within the PBL enhances the equivalent potential temperature of the lower atmosphere, facilitating the accumulation of convective available potential energy (CAPE). Concurrently, the weakening of the inversion at the top of the PBL and the upward moisture transport reduce convective inhibition (CIN), ultimately triggering convection (Dai, 2024; Dai et al., 1999). However, a characteristic lag of 3–6 h exists between peak surface fluxes and the precipitation maximum. This delay is attributed to the initial development of shallow cumulus, which experience lateral and cloud-top entrainment of dry environmental air into the convective plumes. This process dilutes the buoyancy within the core of the rising air and suppresses its vertical growth (Tian and Kuang, 2016). Consequently, a transitional period is required for shallow convection to moisten the lower free atmosphere and further erode CIN. As the transition proceeds and the plumes grow, a smaller fraction of the plume is exposed to the environment and the entrainment rate itself decreases, so the dilution weakens and deep convection can develop (Stirling and Stratton, 2012). Through these complex land-atmosphere interactions, the precipitation maximum over most terrestrial regions is realized in the late afternoon to early evening, between 15:00 and 21:00 local solar time (LST) (Dai et al., 1999; Watters et al., 2021).

However, this surface-heating-driven mechanism is not universally dominant. In regions where large-scale dynamical forcing governs precipitation, diurnal characteristics diverge markedly from the solar-driven pattern. For instance, over the US Great Plains (GP), the nocturnal low-level jet (LLJ) intensifies after sunset due to boundary layer decoupling and transports moist air northward from the Gulf of Mexico. The resulting convergence at the poleward terminus of the jet triggers elevated convection and the propagation of mesoscale convective systems (MCSs), leading to a precipitation peak between midnight and early morning (Higgins et al., 1997; Shapiro et al., 2016). Over the Maritime Continent, daytime convection initiated over land propagates offshore via land breezes and gravity waves, resulting in nocturnal-to-early morning rainfall over coastal waters (Love et al., 2011; Lu et al., 2021; Mori et al., 2004). Furthermore, diurnal cycles differ between convective (CP) and large-scale precipitation (LSP); LSP generally exhibits a weak diurnal amplitude with a dawn peak (approximately 06:00 LST), whereas CP displays a robust diurnal cycle peaking in the late afternoon (16:00–20:00 LST). Since total precipitation (TP) is the aggregate of these components, the timing of the TP peak often precedes that of CP by 1–2 h, reflecting the complex interplay between different precipitation regimes (Dai, 2001; Dai and Trenberth, 2004).

In contrast to terrestrial regions, the oceanic diurnal precipitation maximum typically occurs in the early morning (04:00–08:00 LST), with a diurnal amplitude generally weaker than that over land (Dai, 2001; Dai et al., 2007). Over coastal waters, strong daytime heating over adjacent landmasses drives a land–sea breeze circulation that induces low-level convergence over the ocean during the night, resulting in an early-morning precipitation peak (Deser and Smith, 1998). In the tropical and open oceans where continental influence is minimal, differential longwave radiative cooling between cloudy and clear regions during the night enhances nocturnal destabilization and low-level convergence (Gray and Jacobson, 1977). Furthermore, spatial phase differences in semi-diurnal pressure tides can further intensify this early-morning convergence, playing a crucial role in the diurnal timing of oceanic precipitation (Dai et al., 1999).

Given that the diurnal cycle of precipitation results from the complex coupling of radiation, convection, and boundary layer processes, it serves as a critical benchmark for evaluating the physical fidelity of climate models (Cui et al., 2021; Watters et al., 2021). While climate models can reasonably reproduce precipitation characteristics at seasonal and interannual scales, they consistently struggle to simulate a realistic diurnal cycle (Christopoulos and Schneider, 2021; Dong et al., 2023). For instance, the Coupled Model Intercomparison Project (CMIP) Phase 5 and Phase 6 models typically produce a diurnal peak 3–5 h earlier than observations over land, accompanied by a substantial inter-model spread in amplitude. Over the ocean, although the simulated phase is more comparable to observations, the amplitude remains substantially underestimated, a persistent bias across multiple CMIP generations (Christopoulos and Schneider, 2021; Tang et al., 2021). Many models employ CAPE-based convective triggers that activate convection prematurely in response to boundary layer heating. Furthermore, because these schemes often diagnose shallow and deep convection independently, they fail to represent the gradual transition between these two stages (Freitas et al., 2021; Yang et al., 2021). This leads to the well-known “too frequent, too light” problem in diurnal precipitation and an underestimation of the diurnal cycle amplitude (Dai, 2006; Stephens et al., 2010). Additionally, coarse spatial resolution hinders the proper representation of MCS initiation and propagation, particularly in regions dominated by nocturnal peaks (Song et al., 2024). Reanalysis datasets inherit these systemic problems, as precipitation fields largely rely on model physics with limited sub-daily data assimilation, leading to the same premature peaks found in climate models (Dai and Trenberth, 2004; Watters et al., 2021).

Satellite remote sensing has been widely utilized to address the limitations inherent in climate models and reanalyses. In particular, the Tropical Rainfall Measuring Mission (TRMM) launched in 1997, followed by the Global Precipitation Measurement (GPM) mission, has provided unprecedented global observations of the vertical structure and spatiotemporal variability of precipitation (Kummerow et al., 1998). These high-resolution datasets serve as an essential benchmark for validating model physics and reanalysis performance. With the proliferation of diverse satellite products (e.g., the Integrated Multi-satellite Retrievals for GPM, IMERG; the Climate Prediction Center Morphing Technique, CMORPH; and the Global Satellite Mapping of Precipitation, GSMaP), numerous intercomparisons of the diurnal precipitation cycle have been conducted (Afonso et al., 2020; Janowiak et al., 2005). However, differences in retrieval algorithms and temporal sampling strategies among these products lead to large discrepancies in their captured diurnal characteristics (Pradhan et al., 2025; Watters et al., 2021). Furthermore, because satellite-based precipitation estimations rely on indirect proxies (e.g., radiative emission by liquid at lower levels, radiative scattering by ice at upper levels, and cloud-top temperatures) rather than direct surface measurements, they frequently exhibit a systematic phase delay in the diurnal peak relative to ground-based observations (Guilloteau et al., 2018; Petković et al., 2018; Sun et al., 2018).

As discussed above, both satellite and reanalysis datasets possess inherent biases in their representation of the diurnal precipitation cycle. Although various ground-based studies have attempted to diagnose these issues, several critical gaps remain. First, global-scale diurnal cycle analyses have relied on synoptic present weather codes (WMO ww codes from SYNOP reports) rather than direct precipitation measurements (Dai, 2001; Dai et al., 2007). Moreover, studies using hourly precipitation observations are limited to specific regions or a small number of observation sites (Dai, 2024; Tang et al., 2021), leaving a gap in globally representative, ground-observation-based benchmarks. Second, intercomparisons of satellite products have been confined to comparisons among satellite datasets (Pradhan et al., 2025), and no study simultaneously evaluated the diurnal cycle performance of multiple satellite products against a ground-based observation reference. Third, while phase shifts in individual reanalyses such as the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis Version 5 (ERA5) and the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) have been reported (Dai, 2024; Watters et al., 2021), a comprehensive evaluation remains lacking. Given that convective parameterization has been identified as a primary determinant of diurnal cycle fidelity in global models (Dirmeyer et al., 2012), systematic differences among reanalysis products with distinct convective schemes remain insufficiently explored. This is particularly relevant for the recently improved Japanese Reanalysis for Three Quarters of a Century (JRA-3Q), which features an updated convective parameterization (Kosaka et al., 2024). In particular, the systematic diagnosis of phase shifts through the partitioning of convective and large-scale precipitation remains unexplored across these diverse datasets.

This study addresses these gaps through three objectives. The first is to establish a ground-based benchmark for how precipitation variability is partitioned across temporal frequencies, quantifying the fraction of total variance that resides at sub-daily scales and the proportion of that associated with a diurnal cycle. The second is to evaluate multiple satellite products and reanalyses against this benchmark, extending prior work that has not jointly assessed the diurnal cycle across both product types using a common, dense ground-based reference. The third is to identify the physical origins of the diurnal phase differences among datasets, first by partitioning precipitation into its convective and large-scale components, and then by examining the thermodynamic conditions governing convective initiation in each reanalysis.

To meet these objectives, we use hourly observations from 6824 ground-based stations across Northern Hemisphere (NH) land regions. We evaluate sub-daily precipitation characteristics across five satellite products (IMERG; the TRMM Multi-satellite Precipitation Analysis, TMPA; CMORPH; GSMaP; and the Multi-Source Weighted-Ensemble Precipitation, MSWEP) and four reanalysis datasets (ERA5; JRA-3Q; MERRA-2; and the North American Regional Reanalysis, NARR). Unlike previous studies that focused primarily on the phase and amplitude of the diurnal cycle, this work first quantifies the fraction of TP variability explained by sub-daily fluctuations. We then diagnose the relative contributions of the diurnal (24 h), semi-diurnal (12 h), and higher-order harmonics within the sub-daily band to evaluate the harmonic balance inherent in each product. Furthermore, the study compares the dominant diurnal phase and amplitude across all datasets and investigates the physical drivers of phase shifts, particularly in reanalysis products. Finally, radiosonde observations from the Atmospheric Radiation Measurement (ARM) Southern Great Plains (SGP) site are integrated to explore the boundary layer and convective processes underlying the systematic biases identified in the reanalyses.

The remainder of this paper is organized as follows. Section 2 describes the ground-based observations, satellite products, and reanalysis datasets utilized in this study. Section 3 describes the analytical methods, including sub-daily variance decomposition and harmonic analysis. Section 4 presents the key results, including the decomposition of global precipitation variability across various time scales, an intercomparison of diurnal cycle phase and amplitude, and an analysis of diurnal characteristics for the CP and LSP components. The physical mechanisms are then examined over Northern Hemisphere land under contrasting moisture-supply regimes and at two representative sites, the ARM SGP site and a site in the southeastern United States. Finally, Sect. 5 summarizes the major findings and discusses their broader implications for climate modelling and future observations.

2 Data

2.1 Ground-based observations

2.1.1 Global precipitation observations

Ground-based rain gauge observations provide direct surface measurements that are widely used as the reference for evaluating satellite and reanalysis datasets (e.g., Dai et al., 2007; Tang et al., 2021; Xie et al., 2007). The observational network comprises six distinct data sources (Table 1): (1) a subset of the Global Sub-Daily Rainfall Dataset (GSDR; Lewis et al., 2019), which provides quality-controlled hourly precipitation from approximately 23 000 gauge stations worldwide, (2) the Hadley Centre Integrated Surface Database (HadISD; Dunn et al., 2016), containing sub-daily precipitation observations from synoptic weather stations with comprehensive quality control, hourly precipitation datasets from (3) Germany Deutscher Wetterdienst (DWD), (4) Environment and Climate Change Canada (ECCC), and (5) Korea Meteorological Administration (KMA) national weather stations, and (6) China Meteorological Forcing Dataset version 1.6 (CMFD; He et al., 2020), which provides gridded 3-hourly precipitation at 0.1° resolution over China based on a fusion of 753 China Meteorological Administration station observations with satellite-based data. These precipitation observations have 1- or 3-hourly temporal resolution.

Table 1 Description of ground-based precipitation datasets. The numbers in parentheses in the spatial resolution column indicate the total number of stations used in this study.

The spatial distribution of gauge observations exhibits substantial regional heterogeneity, with station density predominantly concentrated in the NH, particularly over North America, Europe, and East Asia. To ensure sufficient temporal representativeness, only stations with complete 3-hourly records for at least 3 boreal summers (JJA; 92 d × 3 years = 276 d) are retained. The resulting station network is concentrated in NH mid-latitude land regions, reflecting the global heterogeneity of long-term sub-daily gauge availability. We therefore restrict quantitative validation against ground observations to NH land, while presenting global maps in subsequent figures to characterize the spatial structure of inter-product differences. Over ocean and Southern Hemisphere land, where direct validation is not feasible, multi-satellite (MS) and multi-model (MM) means serve as references for relative inter-product comparison rather than absolute accuracy assessment.

To minimize the influence of coastal precipitation regimes, which are governed by land–sea breeze circulations and offshore propagation of convective systems that are not representative of continental diurnal forcing, stations located within 50 km of the coastline are excluded based on their coordinates. The same distance criterion is applied to CMFD grid cells using their centre coordinates. After applying both the coastal exclusion and the minimum record length criteria described above, the final dataset comprises 6071 point-based stations (North America: 3975, Europe: 1675, East Asia: 421) and 753 CMFD station-based gridded observations over China, yielding a combined total of 6824 observation records (Fig. 1). To minimize inconsistencies in spatial representativeness across datasets, the station statistics are spatially averaged onto a 0.5° × 0.625° grid, and CMFD grid cells are remapped to the same grid, corresponding to the coarsest spatial resolution among the gridded products (MERRA-2). The station statistics and the native CMFD cells are aggregated to each target grid cell by simple averaging of all stations or cells whose centre coordinates fall within it. A grid cell is assigned a value whenever at least one valid station or cell is present, with no minimum number imposed, as the statistics are computed at each station prior to spatial aggregation, and partially covered cells at the network edges are retained on the same basis. This procedure is consistent with the spatial harmonization described in Sect. 3.

Figure 1 Locations of ground-based precipitation observations over the Northern Hemisphere. Symbols denote point-based observation stations from each dataset, and the shaded area indicates the spatial coverage of the gridded China Meteorological Forcing Dataset (CMFD) dataset, which is a fusion product combining station observations with reanalysis data. Bracketed numbers indicate the total number of stations or the number of stations to construct the gridded CMFD product. CMFD data courtesy of CMA and collaborators. Boxes delineate the three major regions of station concentration (North America, Europe, and East Asia), with the total number of observations within each region indicated below each box label. Stars mark the two case-study sites examined in Sect. 4.3.2, the Southern Great Plains site (black; 36.6025° N, 97.4925° W) and the southeastern US site (red; 33.5° N, 82° W).

Sub-daily precipitation variance is computed independently at each ground-based station prior to aggregation within each grid cell (hereafter, variance-then-averaging approach). An alternative approach, in which station time series are first averaged at each time step across all contributing stations within a grid cell before variance is computed from the resulting composite time series (hereafter, averaging-then-variance approach), would introduce a systematic dependence on the number of contributing stations (N), since the variance of the mean of random variables is proportional to $\mathrm{1}/N$ of the individual variance. Grid-cell variance estimates from that approach would therefore be biased as a function of station density, an attenuation documented in prior studies (Osborn and Hulme, 1997; Light et al., 2022). The variance-then-averaging approach adopted here follows the “statistic-then-grid” order of operations employed in widely used observational datasets (Donat et al., 2013; Gibson et al., 2019).

Figure 2 confirms this behaviour at a representative grid cell in the Central United States, where station subsets of increasing size N are randomly drawn from the full set of N_max available stations. The variance-then-averaging estimate ($\stackrel{\mathrm{‾}}{{\mathit{\sigma }}^{\mathrm{2}}}$) remains approximately constant across all values of N, whereas the averaging-then-variance estimate (${\mathit{\sigma }}^{\mathrm{2}}\left(\stackrel{\mathrm{‾}}{x}\right)$) decreases monotonically as N increases, at a rate more gradual than the theoretical $\mathrm{1}/N$ scaling because neighbouring stations are positively correlated. The same behaviour is found at grid cells in Germany and East China (Fig. S1 in the Supplement). The spread in precipitation variance among stations within each grid cell indicates how well the grid-cell mean represents station-level values (Fig. S2). The mean absolute deviation remains below 13 % of the grid-cell mean and shows no systematic relationship with the number of stations.

Figure 2 Sensitivity of the grid-cell averaged sub-daily precipitation variance (σ²) to the number of stations (N), evaluated at a representative 0.5° × 0.625° grid cell in the Central United States (33.0° N, 96.9° W). For each N from 1 to N_max, all possible combinations of station subsets are used to compute ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ under two approaches. The $\stackrel{\mathrm{‾}}{{\mathit{\sigma }}^{\mathrm{2}}}$, in which variance is computed at each station individually and subsequently averaged across the selected stations (red line). ${\mathit{\sigma }}^{\mathrm{2}}\left(\stackrel{\mathrm{‾}}{x}\right)$, in which station time series are first spatially averaged and variance is subsequently computed from the resulting mean time series (blue line). The dashed line indicates the theoretical $\mathrm{1}/N$ scaling expected when variance decreases solely as a result of averaging independent samples.

2.1.2 SGP site observations

Observational data from the Atmospheric Radiation Measurement (ARM) Southern Great Plains (SGP) observatory, located near Lamont, Oklahoma (36.6025° N, 97.4925° W, 306 m elevation), are used to investigate the physical processes underlying the diurnal cycle of precipitation. The ARM SGP site provides qualified long-term observations of precipitation, atmospheric state variables, and surface energy fluxes that are essential for understanding the mechanisms of diurnal precipitation variability over the central United States.

Hourly precipitation is obtained from the ARM Best Estimate Atmospheric Measurements (ARMBEATM) product (Chen and Xie, 1994). ARMBEATM provides the best estimate of atmospheric states specifically tailored for evaluating global climate models. The precipitation measurements are obtained from both standard and optical rain gauges at the central facility. High-frequency raw measurements at 30 s or 1 min frequency are temporally averaged to an hourly resolution to ensure consistency with model output. These single-point measurements provide quality-controlled precipitation records suitable for a detailed analysis of the diurnal cycle, although they may not fully capture the spatial variability within gridded model domains.

The ARM Variational Analysis (VARANAL) product (Tang et al., 2019; Zhang et al., 2001) provides large-scale forcing fields derived from ground-based observations constrained by ARM surface and top-of-atmosphere measurements. VARANAL was originally designed to drive single-column models, cloud-resolving models, and large-eddy simulations by providing atmospheric boundary conditions derived from consistent estimates of atmospheric state variables and surface fluxes. The SGP continuous forcing dataset combines two data streams, sgp60varanarapC1.c1 (based on Rapid Refresh model) and sgp60varanarucC1.c1 (based on Rapid Update Cycle model), representing conditions averaged over a domain of approximately 370 km × 300 km centered on the SGP central facility at 1-hourly temporal resolution. In this study, surface sensible heat flux and latent heat flux from VARANAL are utilized to examine their relationship with the diurnal precipitation cycle.

Upper-air atmospheric conditions are characterized using radiosonde observations from the Integrated Global Radiosonde Archive (IGRA) version 2 (IGRA2; Durre et al., 2006, 2016) for station USM00074646 located at Lamont, Oklahoma (36.6025° N, 97.4925° W, 306 m elevation), which provides 6-hourly soundings from 2000 to present. The IGRA2 data are used in the place of upper-air measurements from ARMBEATM to calculate CAPE, CIN, moist static energy (MSE), and atmospheric precipitable water tendency (dPW). All observed variables at the SGP site are utilized for the period 2000–2023, using available days within this period. Surface sensible and latent heat fluxes from VARANAL are used as 3-hourly averages, while precipitation data are used as 3-hourly accumulations.

2.2 Satellite-based datasets

The satellite products and reanalyses evaluated in this study were selected to span the range of design choices relevant to diurnal precipitation characteristics, so that assessing them against direct ground-based observations provides a reference for selecting among them at this time scale. The Integrated Multi-satellite Retrievals for GPM (IMERG) Version 07B Final Run (Huffman et al., 2020, 2023) is a unified algorithm that intercalibrates, merges, and interpolates precipitation estimates from multiple sources including passive microwave (PMW) retrievals from the GPM constellation satellites, microwave-calibrated infrared (IR) satellite estimates, and monthly gauge analyses. IMERG provides global precipitation estimates at 0.1° spatial resolution and 30 min temporal resolution from 1998 to present, using the GPM Core Observatory and TRMM satellites as calibration references.

The algorithm propagates PMW-retrieved precipitation fields forward and backward in time through quasi-Lagrangian interpolation, in which the motion vectors are derived from precipitation (PRECTOT), total precipitable liquid water (TQL), and total column water vapor (TQV) fields of the MERRA-2 analyses (Huffman et al., 2026). A Kalman filter subsequently merges these propagated fields with the estimates derived from the global collection of geostationary IR brightness temperatures (Janowiak et al., 2001) to yield weighted half-hourly precipitation estimates (Joyce and Xie, 2011). The Final Run incorporates monthly gauge adjustments from the Global Precipitation Climatology Centre (GPCC), using the Version 2022 gauge analysis (Schneider et al., 2014, 2022), applying spatially varying ratio multipliers that remain constant throughout each month to correct the bias in satellite-based estimates.

The Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA) research product (3B42) Version 7 (Huffman et al., 2007) provided 3-hourly precipitation estimates at 0.25° × 0.25° spatial resolution between 50° N and 50° S from 1998 to 2019 when its production was discontinued. TMPA combined PMW precipitation estimates from multiple low-Earth-orbit satellites with geostationary IR estimates (Table 2). The individual sensor PMW estimates were intercalibrated by histogram matching to a common reference, the TRMM Combined Instrument, which merges the TRMM Microwave Imager (TMI) and Precipitation Radar (PR), and merged into a single microwave field. Geostationary IR brightness temperatures were then converted to precipitation through spatially varying coefficients calibrated against this merged microwave field, and these IR-based estimates filled the gaps between microwave overpasses. The research product additionally applied a monthly gauge adjustment, scaling the satellite estimates to match the monthly accumulations of the GPCC gauge analysis (Schneider et al., 2014). Because this adjustment applied a single monthly ratio to every 3-hourly field within the month, it constrained the long-term accumulation while preserving the relative sub-daily structure. TMPA is retained despite its discontinued production because, as the predecessor to IMERG, it fills temporal gaps using calibrated infrared estimates rather than propagating microwave retrievals in time. This contrast with IMERG, CMORPH, and GSMaP allows the effect of temporal propagation on sub-daily variability to be assessed across algorithm generations.

The Climate Prediction Center (CPC) Morphing Technique (CMORPH; Joyce et al., 2004) produces global precipitation estimates by integrating PMW precipitation retrievals from multiple low earth orbit satellites (TMI on TRMM, AMSR on AQUA, SSM/I and SSMIS on DMSP F-13–F-18, AMSU and Microwave Humidity Sounder (MHS) on NOAA-15–19 and MetOp-A/B, and MWRI on FY-3B). Unlike other satellite products that directly combine PMW and IR estimates, CMORPH uses cloud motion vectors derived from consecutive geostationary IR brightness temperature images to realize the propagation of PMW retrievals through time and space. The morphing technique employs a time-weighting scheme that blends forward- and backward-propagated PMW estimates, with weights inversely proportional to the temporal distance from the observation time. This approach yields precipitation estimates at 30 min intervals, effectively bridging the gaps between infrequent PMW overpasses. This morphing approach provided the conceptual basis for the propagation scheme later adopted in IMERG (Joyce and Xie, 2011), although the two algorithms have since diverged, with IMERG now deriving its motion vectors from reanalysis fields rather than from geostationary IR (Tan et al., 2019).

This study uses the bias-corrected CMORPH Climate Data Record (CDR) Version 1.0 (Xie et al., 2017), which applies systematic bias correction to the raw satellite estimates. Over land, bias correction is performed through daily-based probability density function (PDF) matching against the CPC Unified Daily Gauge Analysis, which adjusts the statistical distribution of satellite-estimated precipitation rates to match ground-based observations. Over ocean, the Global Precipitation Climatology Project (GPCP) Pentad (5 d) Product is used to carry out the PDF matching technique in two steps: first using climatological correction based on a 15-year dataset, then applying real-time refinement using a 31 d moving window (Xie et al., 2017). The adjustment coefficients, calculated as the ratio between the mean reference field (CPC gauge or GPCP pentad) and the mean raw CMORPH, are downscaled from daily 0.25° resolution to the native 30 min 8 km resolution. The dataset has a global coverage from 60° S to 60° N spanning from January 1998 to the present.

The Global Satellite Mapping of Precipitation (GSMaP) version 8 gauge-calibrated product (GSMaP_Gauge) provides global hourly precipitation estimates at 0.1° × 0.1° spatial resolution between 60° N and 60° S (Kubota et al., 2020). GSMaP retrieves precipitation using PMW radiometer data from the GPM constellation and integrates them with IR observations from geostationary satellites via a Kalman filtering algorithm, which propagates PMW estimates both forward and backward in time (Ushio et al., 2009). The gauge-calibrated version further refines these satellite-based estimates by incorporating daily rain gauge data from the National Oceanic and Atmospheric Administration (NOAA) CPC Unified Daily Gauge Analysis to enhance accuracy over land (Kubota et al., 2022; Lu and Yong, 2020; Mega et al., 2019). Version 8, released in December 2021, introduces substantial improvements including updated PMW retrieval databases, enhanced orographic rainfall estimation, and a normalization module to homogenize retrievals among multiple microwave sensors (Gao et al., 2024; Tashima et al., 2021; Yamamoto and Kubota, 2022). Additionally, the retrieval domain has been extended toward the poles during PMW overpasses.

The Multi-Source Weighted-Ensemble Precipitation (MSWEP) Version 2.8 (Beck et al., 2019) is a merged dataset that integrates gauge observations, satellite products, and reanalysis outputs using a weighted-ensemble methodology. At a 3 h time step and 0.1° spatial resolution, the non-gauge component combines reanalysis output (ERA5) with satellite-based estimates, for which IMERG is used after 2000 and infrared-based estimates derived from GridSat-B1 brightness temperatures are used for the pre-2000 period. The GridSat-based precipitation is obtained by cumulative distribution function matching the infrared brightness temperatures to CMORPH within air-temperature zones, rather than taken as a direct precipitation product. Weights are assigned based on station density for gauge-based components and comparative skill against nearby gauges for satellite and reanalysis components. All possible combinations of datasets are merged through weighted averaging, with cumulative weights calculated to represent the total information content. The merged satellite–reanalysis estimates of each dataset combination are then cumulative distribution function matched at the 3-hourly time scale to a reference distribution based on the Multi-Source Weather (MSWX) product (Beck et al., 2022) to mitigate spurious light rainfall and the underestimation of heavy precipitation.

MSWEP incorporates systematic daily gauge corrections utilizing multiple sources including the Global Historical Climatology Network–Daily (GHCN-D) and the Global Summary of the Day (GSOD) through a multiplicative approach that preserves sub-daily distributions. Notably, the dataset accounts for regional variations in gauge reporting times (ranging globally from 00Z to 12Z), thereby minimizing temporal mismatches between gridded estimates and gauge observations. Monthly climatology corrections are applied using the Climatologies at High Resolution for the Earth's Land Surface Areas (CHELSA; Karger et al., 2017), replacing the WorldClim-based corrections used in earlier versions. MSWEP further corrects systematic terrestrial biases (e.g., gauge undercatch and orographic effects) by inferring catchment-average precipitation from river discharge observations at 13 762 catchments globally, a discharge-based scheme that distinguishes it from the gauge-based undercatch corrections applied in other products (Beck et al., 2019). The dataset provides global coverage at 0.1° and 3-hourly resolution from 1979 onward.

Table 2 Description of satellite-based precipitation datasets. Primary data sources include passive microwave (PMW) and infrared (IR) sensors. Temporal resolution of bias correction is indicated in parentheses (D: daily; 5D: 5 d pentad; M: monthly).

2.3 Reanalysis datasets

The European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis Version 5 (ERA5; Hersbach et al., 2020) provides global hourly atmospheric estimates from 1950 to present. It employs the Integrated Forecasting System (IFS) Cycle 41r2 with a hybrid incremental four-dimensional variational (4D-Var) data assimilation scheme, assimilating observations with 12 h assimilation windows (09:00–21:00 UTC and 21:00–09:00 UTC). The atmospheric model has 0.25° × 0.25° horizontal resolution with 137 hybrid sigma-pressure vertical levels extending from the surface to 0.01 hPa. ERA5 reanalysis incorporates the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (H-TESSEL) to simulate land surface variables, utilizing a weakly coupled land data assimilation system to adjust soil moisture and snow depth with screen-level temperature and humidity observations (Balsamo et al., 2009). Notably, ERA5 precipitation is derived from 18 h forecasts initialized twice daily (06:00 and 18:00 UTC), with hourly outputs accumulated from the first 12 h (1–12 h lead time) (Lavers et al., 2022; Soci et al., 2024). Each hourly value corresponds to a forecast lead time of 1–12 h, depending on its temporal proximity to the initialization. Crucially, precipitation observations are not directly assimilated; instead, precipitation is generated as a diagnostic output from the model physics.

The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), is NASA's global atmospheric reanalysis produced by the Global Modeling and Assimilation Office (GMAO), providing data from 1980 to the present at 0.5° × 0.625° horizontal resolution with hourly temporal resolution (Gelaro et al., 2017). The system utilizes the Goddard Earth Observing System Model, Version 5 (GEOS-5; Molod et al., 2015) with a 6-hourly assimilation window based on the Gridpoint Statistical Interpolation scheme (Wu et al., 2002) and Incremental Analysis Update (IAU) methodology (Bloom et al., 1996). The assimilation system incorporates extensive conventional and satellite-based atmospheric observations, including modern microwave sounders such as the Advanced Microwave Sounding Unit (AMSU-A), Advanced Technology Microwave Sounder (ATMS), and MHS that provide atmospheric temperature and humidity profiles, and hyperspectral infrared radiance instruments such as the Atmospheric Infrared Sounder (AIRS), Infrared Atmospheric Sounding Interferometer (IASI), and Cross-track Infrared Sounder (CrIS) that provide temperature and moisture sounding information. Precipitation fields are generated as a diagnostic output from model physics, without direct correction or assimilation of precipitation observations (Table 3).

The Japanese Reanalysis for Three Quarters of a Century (JRA-3Q) is the third global atmospheric reanalysis produced by the Japan Meteorological Agency (JMA), covering September 1947 to present at TL479 (∼ 0.375°) horizontal resolution (Kosaka et al., 2024). JRA-3Q builds upon its predecessor Japanese 55-year Reanalysis (JRA-55; Kobayashi et al., 2015), by incorporating three major updates: (1) an updated forecast model and data assimilation system based on the JMA operational numerical weather prediction system since December 2018, (2) improved boundary conditions, including updated sea surface temperature and sea ice datasets, and (3) the assimilation of additional observations including Global Navigation Satellite System (GNSS) zenith total delay (ZTD) for atmospheric water vapor, GNSS radio occultation (up to 60 km altitude), airport surface pressures, and hyperspectral infrared sounder radiances (e.g., AIRS, IASI, CrIS) for temperature and humidity profiles. The reanalysis employs a 4D-Var data assimilation system using the JMA Global Spectral Model (GSM), with assimilation cycles every 6 h (00:00, 06:00, 12:00, and 18:00 UTC).

The North American Regional Reanalysis (NARR), produced by NOAA's National Centers for Environmental Prediction (NCEP), provides regional atmospheric data over North America from 1979 to present at ∼ 32 km horizontal resolution and 3-hourly temporal resolution (Mesinger et al., 2006). NARR employs the NCEP Eta Model with a three-dimensional variational (3D-Var) data assimilation system, cycling every 3 h where the 3 h forecast from the previous cycle serves as the first guess. Unlike global reanalyses, NARR assimilates precipitation observations via a latent heating nudging technique (Lin et al., 1999). This method converts observed precipitation into vertical profiles of latent heating, which are then nudged into the model's temperature tendency equation during the forecast step. Over the continental United States, the system assimilates the Hourly Precipitation Dataset (HPD) and Parameter-elevation Regressions on Independent Slopes Model (PRISM; Daly et al., 1994) daily precipitation incorporating diurnal patterns from the Global Reanalysis 2 (GR2) forecasts when the HPD data is unavailable. For Canada and Mexico, GR2-based diurnal patterns are similarly applied to daily national gauge data. Over oceanic regions (south of 27.5° N), Climate Prediction Center Merged Analysis of Precipitation (CMAP) pentad data (Xie and Arkin, 1997) are assimilated; however, north of 42.5° N, precipitation assimilation is omitted due to uncertainties in the CMAP product in high latitudes.

Table 3 Description of reanalysis precipitation datasets (DA: data assimilation).

3 Methods

To ensure consistent intercomparison across all precipitation datasets, we perform both spatial and temporal harmonization. All datasets, regardless of their original format (gridded products, point observations, or reanalysis outputs), are remapped to a common spatial resolution of 0.5° × 0.625° (matching the MERRA-2 grid) using bilinear interpolation. This method efficiently transforms heterogeneous inputs onto a unified grid while preserving spatial continuity, thereby minimizing scale-dependent representativeness errors (Gervais et al., 2014; Tustison et al., 2001).

All precipitation products are standardized to 3-hourly accumulated totals (mm 3 h⁻¹). Datasets provided as precipitation rates (e.g., mm s⁻¹, kg m⁻² s⁻¹) are converted to mm h⁻¹ and accumulated over 3 h windows, while those originally supplied as 3-hourly accumulations are used directly. This harmonization protocol ensures that discrepancies between products reflect intrinsic variations in retrieval algorithms and data sources rather than inconsistencies in resolution or units (Pradhan et al., 2025). After this data pre-processing, diurnal cycle analyses are restricted to rainy days (daily precipitation ≥ 0.275 mm d⁻¹) to exclude unreliable signals from weak precipitation signals (Watters et al., 2021).

This study focuses exclusively on the boreal summer season (June–August; JJA), when diurnal precipitation variability is strongest over Northern Hemisphere land. Analysis periods vary due to data availability: CMORPH and GSMaP (1998–2023), TMPA (1998–2019), IMERG (1998–2023) and MSWEP (2000–2023), and all reanalysis products (2000–2023). Previous studies suggest that a stable diurnal cycle can be derived from just a few years of data (Covey et al., 2016; Dai et al., 2009; Tang et al., 2021), supporting the robustness of cross-period comparisons despite differences in temporal coverage among datasets. To further assess the sensitivity of the results to the choice of analysis period, we conducted sensitivity tests by restricting selected datasets to non-overlapping sub-periods (e.g., CMORPH: 1998–2007 and 2013–2022) and comparing the resulting diurnal phase estimates against those derived from the full analysis period. The mean phase difference between sub-periods and the full-period analysis is less than 1 h across all tested datasets, confirming that the diurnal cycle characteristics are not substantially sensitive to the choice of analysis period (not shown).

3.1 Sub-daily variance decomposition

To quantify the relative contribution of sub-daily (period < 24 h) variations to the TP variability, a 30th-order Butterworth high-pass filter with a cutoff period of 26 h (Oppenheim and Schafer, 2010) is applied, separating sub-daily frequency components from relatively lower-frequency components (period ≥ 24 h, including daily mean), where their respective variances are denoted by ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{Low}}^{\mathrm{2}}$. The mathematical formulation of the sub-daily variance ratio (rat_subD) is followed as:

$$\begin{array}{}\text{(1)}& {\mathrm{rat}}_{\mathrm{subD}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}}}\times \mathrm{100}\mathit{\%}\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ is the variance of the unfiltered time series.

Following Dong et al. (2023) and Covey et al. (2018), ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$is further decomposed into the regular sub-daily variance (${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$), associated with the climatological mean diurnal cycle, and the irregular sub-daily variance (${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$), representing sub-daily fluctuations that deviate from the climatological mean diurnal cycle on individual days:

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}\end{array}$$

The regular component is defined as the variance of the climatological mean diurnal cycle $P_i^{\mathrm{reg}}$, obtained by averaging the high-pass filtered precipitation time series over all days at each time of day i:

$$\begin{array}{}\text{(3)}& {\displaystyle }{P}_i^{\mathrm{reg}}={\displaystyle \frac{\mathrm{1}}{D}}{\sum }_{n=\mathrm{1}}^D{P}_{i,n},i=\mathrm{1},\mathrm{\dots },N\text{(4)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\left(P_i^{\mathrm{reg}}-\stackrel{\mathrm{‾}}{P^{\mathrm{reg}}}\right)}^{\mathrm{2}}\end{array}$$

where N is the number of sub-daily time steps per day (N=8 for 3-hourly temporal resolution), D is the total number of days, and $\stackrel{\mathrm{‾}}{P^{\mathrm{reg}}}$ is the daily mean of $P_i^{\mathrm{reg}}$. The irregular component is obtained as the residual:

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}\end{array}$$

The regular component is further decomposed into Fourier harmonics via Fast Fourier Transform. Since $P_i^{\mathrm{reg}}$is a periodic sequence with period N, its variance decomposes as:

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2}\mathrm{nd}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{higher}}^{\mathrm{2}}\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{2}\mathrm{nd}}^{\mathrm{2}}$ correspond to the diurnal (24 h), semidiurnal (12 h) cycles, respectively, and ${\mathit{\sigma }}_{\mathrm{higher}}^{\mathrm{2}}$ is the residual variance from higher-order harmonics (period ≤ 8 h).

3.2 Characterization of regular sub-daily components

UTC-based time coordinates are converted to Local Solar Time (LST; t_LST) to explicitly account for the solar forcing of diurnal precipitation (LST = UTC $+\mathit{\lambda }/\mathrm{15}\mathit{°}$, where λ is the longitude in degrees). As 3-hourly time series require an adjustment to align with the discrete time steps, we first calculate $\mathit{\lambda }/\mathrm{45}\mathit{°}$, round the result to the nearest integer, and subsequently multiply by 3. This conversion ensures that the diurnal cycle is analysed relative to local solar noon rather than a fixed UTC reference. Such alignment is critical for identifying physically meaningful precipitation phases driven by solar heating.

The sub-daily precipitation is characterized through harmonic analysis, where the precipitation time series are decomposed into Fourier components at each grid point (Cui et al., 2021; Dai, 2001):

$$\begin{array}{}\text{(7)}& P\left(t_{\mathrm{LST}}\right)=\stackrel{\mathrm{‾}}{P}+{\sum }_{k=\mathrm{1}}^{N/\mathrm{2}}{A}_k\sin \left(k\mathit{\omega }{t}_{\mathrm{LST}}+{\mathit{\phi }}_k\right)\end{array}$$

where $\stackrel{\mathrm{‾}}{P}$ is the daily-mean precipitation, A_k and φ_k are the amplitude and phase of the k-th harmonic, N is the number of time samples per day (N=8 for 3-hourly data), and $\mathit{\omega }=\mathrm{2}\mathit{\pi }/\mathrm{24}$. For instance, the 1st harmonic (k=1, 24 h period) represents the primary diurnal signal, while the 2nd harmonic (k=2, 12 h period) captures semi-diurnal variability.

To characterize the relative contribution of individual harmonic components within ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$, Eq. (6), the fractional contribution of each harmonic component to ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ is defined as:

$$\begin{array}{}\text{(8)}& {\mathrm{rat}}_k={\displaystyle \frac{{\mathit{\sigma }}_k^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}}}\times \mathrm{100}\mathit{\%}\end{array}$$

We quantify the relative dominance of variance at specific frequency bands within the regular sub-daily variability (e.g., 1st (k=1), 2nd (k = 2), and high-order (k≥3) harmonics), thereby enabling an assessment of the fidelity of various datasets in capturing the observed spectral balance.

To assess the adequacy of the 1st harmonic component in representing the observed diurnal cycle, we use the rat₁. Regions where rat₁ ≥ 50 % indicate that the 24 h periodic component accounts for at least half of ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$. The full hierarchy of the variance decomposition described in Sect. 3.1 and Sect. 3.2 is summarized in Fig. 3.

Figure 3 Schematic of the sub-daily precipitation variance decomposition, showing the method applied at each tier and the corresponding equations in Sect. 3.1 and 3.2. Box widths are uniform and do not indicate the relative magnitude of each component.

3.3 Evaluation of diurnal phase and amplitude

Gridded satellite-based and modelled datasets are evaluated against ground-based observations with respect to the diurnal characteristics of precipitation phase and amplitude. For phase evaluation, we account for the circular nature of the diurnal cycle by calculating the circular mean bias and circular correlation coefficient. The circular mean bias is computed by converting precipitation peak phase into angular values, determining the circular mean difference between the datasets and observation, and subsequently converting the angular difference back to hours. This approach properly resolves periodicity issues, for instance, recognizing that 23:00 LST and 01:00 LST are only 2 h apart rather than 22 h. Furthermore, we assess the spatial representation of diurnal phase in the datasets using the circular correlation coefficient (Fisher and Lee, 1994). Its mathematical formulation is followed as:

$$\begin{array}{}\text{(9)}& r_c={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\sin \left({\mathit{\theta }}_i-\stackrel{\mathrm{‾}}{\mathit{\theta }}\right)\sin \left({\mathit{\varphi }}_i-\stackrel{\mathrm{‾}}{\mathit{\varphi }}\right)}{\sqrt{{\sum }_{i=\mathrm{1}}^n{\sin }^{\mathrm{2}}\left({\mathit{\theta }}_i-\stackrel{\mathrm{‾}}{\mathit{\theta }}\right){\sum }_{i=\mathrm{1}}^n{\sin }^{\mathrm{2}}\left({\mathit{\varphi }}_i-\stackrel{\mathrm{‾}}{\mathit{\varphi }}\right)}}}\end{array}$$

where θ and ϕ represent the diurnal phase, expressed by the azimuthal angle (in radians), from the dataset and observations at each grid cell (i) and the overbar represents the circular mean. The circular correlation coefficient ranges from −1 to +1, where values close to +1 indicate robust agreement in the spatial patterns of the diurnal phase. Statistical significance is assessed via the Fisher-Lee test, which approximates the distribution of the circular correlation under the null hypothesis of no association.

For amplitude evaluation, conventional linear statistics are employed, given that amplitude values are non-circular continuous variables. The mean bias is calculated as the arithmetic mean of differences between the dataset and observed amplitudes across all stations, while spatial correlation is quantified using Spearman's rank correlation coefficient. Spearman's method is preferred over Pearson correlation due to its robustness to outliers and independence from linear assumptions, making it particularly appropriate given the high variability of precipitation amplitudes across different climate regimes.

3.4 Atmospheric water budget framework

To elucidate the physical processes underlying the occurrence of diurnal precipitation, the diurnal variations of atmospheric variables controlling moisture availability and convective triggering are examined. CAPE and CIN are computed with MetPy (May et al., 2022) using a surface-based parcel on the pressure levels common to all four reanalyses, namely 1000, 925, 850, 700, 500, and 400 hPa, with the parcel lifted from the lowest level available at each grid point. Following the atmospheric water budget framework (Song and Wei, 2021), the vertically integrated moisture flux convergence (MFC) is derived as a residual from the water budget equation:

$$\begin{array}{}\text{(10)}& \mathrm{MFC}=P-E+{\mathrm{dPW}}_t\end{array}$$

where E is evapotranspiration, P is precipitation, and dPW_t is the tendency of atmospheric precipitable water (PW), estimated as ${\mathrm{dPW}}_t={\mathrm{PW}}_{t-\mathrm{1}}-{\mathrm{PW}}_t$. For the observational reference at the ARM SGP site, E is derived from the latent heat flux of the VARANAL product, P from ARMBEATM, and PW from IGRA2 radiosonde observations for each reanalysis product, the corresponding variables from each dataset are used.

MSE is also examined as an integrated measure of the thermodynamic state of the atmosphere. Notably, MSE is considered approximately conserved during adiabatic vertical motions and consists of sensible and latent heat, and geopotential energy, followed as:

$$\begin{array}{}\text{(11)}& \mathrm{MSE}=C_{\mathrm{p}}T+L_{\mathrm{v}}q+gz\end{array}$$

where C_p is specific heat at constant pressure, T is temperature, g is gravitational acceleration, z is geopotential height, L_v is the latent heat of vaporization, and q is specific humidity. MSE is computed over the range 925–400 hPa, effectively representing the thermodynamic conditions of the lower-to-middle troposphere essential for convective triggering.

4 Results

4.1 Precipitation variability in sub-daily time scales

4.1.1 Partitioning of the total variance

Total variance of precipitation at observational sites over NH land during boreal summer is 0.87 mm² h⁻² (Fig. 4a), whereas most precipitation datasets underestimate this variance, and the underestimation is more apparent in the reanalysis products. In contrast, TMPA overestimates the variance. It retains PMW estimates where available and fills the remaining grid boxes with IR estimates calibrated against the PMW estimates (Huffman et al., 2007). The data source at a given grid box can therefore switch between PMW and IR across successive 3-hourly time steps, potentially introducing source-dependent differences in the statistical properties of the resulting time series. This source switching, together with variability in the TRMM Combined Instrument estimates used both as input and as the calibration reference, as well as in the PMW retrievals from the Goddard Profiling Algorithm and the IR estimates, likely contributes to the excess variance. The underestimation of total variance is primarily attributed to the underestimated ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ despite the underestimation of ${\mathit{\sigma }}_{\mathrm{Low}}^{\mathrm{2}}$. For instance, while ${\mathit{\sigma }}_{\mathrm{Low}}^{\mathrm{2}}$ from satellite-based datasets is comparable to observations, they generally underestimate ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$, with the notable exception of TMPA. For reanalysis products, ERA5 reproduces only 64 % of the observed ${\mathit{\sigma }}_{\mathrm{Low}}^{\mathrm{2}}$ and a mere 28 % of the observed ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$. In contrast, MERRA-2 exhibits the most substantial underestimation, capturing only 46 % of the observed ${\mathit{\sigma }}_{\mathrm{Low}}^{\mathrm{2}}$ and 13 % of the ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ relative to observations.

Figure 4 (a) Sub-daily (pink) and lower-frequency (orange) variance components of precipitation at observation locations. Sub-daily frequency variance is obtained via high-pass Butterworth filter. Numbers inside bars indicate variance magnitude mm² h⁻², and percentage above bars represent the sub-daily variance ratio (rat_subD). Differences from the ground-based observations, for (b, d) the Multi-Satellite (MS; IMERG, TMPA, CMORPH, GSMaP, MSWEP) and (c, e) the Multi-Model (MM; ERA5, JRA-3Q, MERRA-2) mean, in terms of (b, c) total variance (${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$) and (d, e) rat_subD, where the spatial averaged value is noted in each panel below.

Further decomposing ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ into regular (${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$) and irregular (${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$) components provides additional insight into the source of this underestimation. In observations, ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$ accounts for approximately 98 % of ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$, indicating that the repeatable diurnal signal is overwhelmed by day-to-day variability in sub-daily precipitation. While convective intermittency and irregular occurrence of MCSs are likely dominant contributors over land (Covey et al., 2018), other sources (e.g., synoptic-scale forcing variability) also project onto ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$. ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ also exhibits considerable deviation from observations across all datasets. Notably, despite the more severe underestimation of ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$ in reanalysis products relative to satellite-based datasets, the relative contribution of ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ within ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ is overestimated in reanalysis products (5 %–8 %) than in satellite-based datasets (2 %–3 %). The disparity in underestimation between ${\mathit{\sigma }}_{\mathrm{Low}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$, and particularly the underrepresentation of ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$, is primarily attributed to two factors. First, satellite-based datasets, where gauge-based observations are typically assimilated to correct relatively low-frequency variations (e.g., daily, pentad, and monthly), leave sub-daily scales less constrained. Second, for reanalysis products, the underestimation of ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ is further exacerbated by the inherent limitations of parameterization schemes in capturing precipitation physics at relatively coarse spatial resolutions (Herold et al., 2016; Skok et al., 2016). These constraints lead to an excessive frequency of light precipitation, which hinders the simulation of intense rainfall during sub-daily time scales (Chen et al., 1996; Dai, 2006; Stephens et al., 2010).

The spatially aggregated statistics in Fig. 4a do not indicate where these discrepancies arise, so the Multi-Satellite (MS) and Multi-Model (MM) means are differenced against the ground-based observations over the regions of dense gauge coverage (Figs. 4b–e). Both product groups underestimate ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ over most of Contiguous United States (CONUS), and the deficit is largest over the central and eastern plains, where it exceeds 1.5 mm² h⁻² in MM (Fig. 4b, c). Warm-season precipitation in this region is dominated by organized convective systems that initiate over the Rockies and propagate eastward overnight (Houze, 2004), so an inability to reproduce the intensity of these systems maps directly onto a deficit in total variance. East Asia separates the two product groups more clearly than CONUS: MS remains broadly consistent with the observations across China, whereas MM still exhibits underestimation over southeastern China. The spatial root-mean-square error (RMSE) summarizes this contrast, at 0.86 mm² h⁻² for MS against 1.19 mm² h⁻² for MM. For the rat_subD, both MS and MM underestimate the observations across nearly the entire valid domain. Spatial RMSE values are 11.46 % for MS and 22.13 % for MM (Fig. 4d, e), with MM showing both a larger magnitude and a more spatially extensive deficit. The ratio of the MM to MS spatial RMSE rises from 1.4 for ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ to 2.0 for rat_subD, indicating that reanalysis products depart from the observations more in the frequency partitioning of the variance than in its total magnitude.

Figure 5 (a) Violin plots of the rat_subD distribution for each dataset over the globe (black, left side), land (red, right side), and ocean (blue, right side), where the thick solid line and circle indicate the median and mean, respectively, and dashed lines denote the 25th and 75th percentiles. Spatial distribution of rat_subD from (b) MS and (c) MM, where crossed areas indicate rat_subD is below 30 %. The spatial correlation of both products against ground-based observations is denoted in the lower left of each panel. The zonal averaged distribution over the land (red) and ocean (blue) is displayed on the right of each panel and shading represents ±1 standard deviation.

The observed rat_subD, quantifying the relative contribution of sub-daily precipitation variability against its total variance, exhibits that 65.6 % of total variance is explained by both systematic diurnal and irregular high-frequency fluctuations. This result serves as a benchmark for evaluating whether datasets accurately reproduce the temporal partitioning of precipitation variability. All satellite-based and reanalysis datasets underestimate the observed rat_subD, with modeled products consistently falling below 50 % of the observed value (ERA5: 45.2 %, JRA-3Q: 40.4 %, MERRA-2: 36.0 %). Regarding the land-ocean contrast, rat_subD over land exceeds that over the ocean in most datasets (Fig. 5a). Most land regions exhibit relatively strong sub-daily precipitation variability relative to the ocean, reflecting greater intensity and frequency of convective precipitation events over land (Li et al., 2025). The lower tail of the negatively skewed distribution corresponds to areas where large-scale dynamics or orographic effects suppress typical diurnal signals. Conversely, oceanic distributions are positively skewed, indicating that while most regions exhibit low rat_subD values reflecting the relatively larger contribution of day-to-day precipitation variability to total variance over the ocean (Covey et al., 2018), the upper tail corresponds to specific areas where reduced low-frequency variability yields disproportionately high sub-daily fractions. GSMaP stands out as a notable exception, exhibiting a reversed land–ocean contrast in rat_subD with a lower median over land than over ocean. This may stem from a combination of two land-specific characteristics of the GSMaP retrieval and calibration system. First, GSMaP is known to exhibit a systematic intensity-dependent bias over land, which overestimates light precipitation frequency while underestimating heavy precipitation intensity (Chen et al., 2019; Ning et al., 2017; Pradhan et al., 2025). It suppresses the physical amplitude of both the regular diurnal cycle and the episodic high-frequency fluctuations relative to other satellite products. Second, the gauge calibration of GSMaP is applied exclusively over land, where the hourly rates are optimized so that their 24 h accumulation matches the daily CPC gauge totals (Mega et al., 2019). Because this optimal-estimation scheme adjusts each hour unevenly rather than applying a single daily multiplier, it can reshape the relative magnitudes of the sub-daily fluctuations over land, while the ocean, lacking gauge constraints, retains its original sub-daily structure. In contrast, IMERG adjusts its half-hourly estimates so that they sum to the monthly satellite-gauge combination, applying a single ratio multiplier that is fixed throughout the month (though spatially varying) to every half-hour within that month (Huffman et al., 2026). Because the same factor scales all half-hours equally, the relative sub-daily structure is preserved. Similarly, CMORPH adjusts distributional shape through probability density function matching without constraining hourly magnitudes (Xie et al., 2017).

The rat_subD spatial patterns from MS and MM reveal distinct regional characteristics (Fig. 5b, c). Both product suites represent the pronounced geographic pattern of strong sub-daily variability in land regions where rat_subD reaches approximately 50 %. The reanalysis products exhibit globally lower rat_subD compared to satellite-based datasets, particularly pronounced over the ocean. In the zonally averaged distribution, satellite-based products show relatively uniform rat_subD meridionally both land and ocean areas. In contrast, the modelled distribution displays a pronounced meridional gradient, characterized by distinct disparities between land and ocean, particularly in the tropics. In terms of spatial agreement with ground-based observations over NH land, the MS demonstrates superior performance over the MM, with a spatial correlation of 0.78 compared with 0.69.

4.1.2 Harmonic composition of the regular component

Although ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$ dominates ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ across all datasets (Fig. 4a), the regular component ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ is of particular diagnostic value as it encodes the physically interpretable periodic structure of sub-daily precipitation variability ­­ including the diurnal and semi-diurnal cycles that are directly linked to surface thermal forcing, boundary layer dynamics, and atmospheric tides. The absolute magnitude of each harmonic component of ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ reveals more pronounced differences than the composition ratios (Fig. 6a), with CMORPH (0.0104 mm² h⁻²) showing the closest agreement with observations (0.0094 mm² h⁻²) among satellite-based products. Among the reanalysis products, all show smaller absolute values of ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ than the observations, with MERRA-2 exhibiting the lowest value (0.0046 mm² h⁻²). ERA5 shows a larger contribution of ${\mathit{\sigma }}_{\mathrm{2}\mathrm{nd}}^{\mathrm{2}}$ relative to other datasets.

Figure 6 (a) Same as Fig. 4a, but for regular sub-daily precipitation variance (${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$) decomposed into the 1st (pink), 2nd (orange), and higher-order (blue) harmonic components. (b–e) Same as Fig. 4b–e, but for (b, c) the 1st harmonic variance (${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$) and (d, e) its fractional contribution to ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ (rat₁).

The regional structure of the ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ difference partly follows that of ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$. Both MS and MM underestimate ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ from the GP to the southern CONUS, with MM showing the larger deficit over a spatially more extensive area (spatial RMSE of 0.015 and 0.017 mm² h⁻², respectively), indicating that the nocturnal convective systems noted above are underestimated in their diurnal amplitude and not only in their total variance. Over China, MS overestimates ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ over southern China, while MM overestimates it more extensively across central and southeastern China, even though ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ is underestimated there (Fig. 4c). The 1st harmonic is therefore not uniformly underestimated across all gridded products. For the rat₁, both MS and MM overestimate it over most of the validation domain, with spatial RMSE values of 17.20 % for MS and 26.61 % for MM, and with the largest values over interior and northern China. An overestimated rat₁ indicates that ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ is carried almost entirely by the diurnal harmonic, so the second and higher-order components remain weaker than observed and the diurnal cycle is represented as smoother than it actually is. The GP is the exception, where MM underestimates rat₁ over the same area in which ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ is most underestimated (Fig. 6c, e).

Figure 7 (a–c) Spatial distributions of ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ and (d–f) ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ from (a, d) MS, (b, e) MM, and (c, f) their difference (MM − MS). The zonal mean distribution over the land (red) and ocean (blue) is displayed on the right of each panel, where the shading represents ±1 standard deviation.

The magnitude of ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ in MS (Fig. 7a) peaks around 30° N over continental regions, with prominent values over central North America, India, southern China, and interior Africa, whereas over the ocean it is generally lower. In contrast to rat_subD (Fig. 5b, c), MM shows generally larger ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ than the MS (Fig. 7b), with values exceeding 10 % near the equatorial land and decreasing poleward before increasing again around 30° N. This arises from the severe underrepresentation of ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$ (cf. Fig. 4a), which reduces ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ disproportionately and thereby inflates the relative contribution of ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$.

The spatial distributions of ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ (Fig. 7d, e) show maximum values over 0–30° N continental regions in both product suites, with particularly prominent signals over southern China, the Amazon, the GP, and the Maritime Continent. The MM (Fig. 7e) consistently underrepresents ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ over both land and ocean, reflecting the characteristic of reanalysis products wherein the relative contribution of ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ within ${\mathit{\sigma }}_{\mathrm{subD}}^{\mathrm{2}}$ is overestimated while its absolute magnitude is underrepresented. This underrepresentation of ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ is not confined to reanalysis products; satellite-based datasets such as GSMaP and MSWEP also exhibit notably lower absolute values relative to observations, suggesting that the physical intensity of the regular diurnal cycle component remains challenging across multiple product types.

On the other hand, the difference between MS and MM in both ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ (Fig. 7c) and ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ (Fig. 7f) reveals that the origin of this overestimation varies by region. For instance, over the western coast of Chile and southern Africa, the MM overestimates ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ despite negligible differences in ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ itself, indicating that the inflated ratio stems from a reduction in ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ driven by the underrepresentation of ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$ rather than from errors in the diurnal cycle amplitude. Over the eastern coast of India, the Maritime Continent, and tropical continental interiors, the MM simultaneously overestimates ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ while underestimating ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$, implying that the far more severe underrepresentation of ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$ reduces ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ to a greater extent than the reduction in ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$. In contrast, over the GP and northwestern India, both ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$ are underestimated by the MM, suggesting that in these regions the insufficient physical intensity of the diurnal cycle itself is the dominant source of error, rather than the underrepresentation of ${\mathit{\sigma }}_{\mathrm{irreg}}^{\mathrm{2}}$.

Figure 8 (a) Ternary diagram of harmonic amplitude ratios for 1st (bottom axis), 2nd (right axis), and higher-order (left axis) harmonic components of regular sub-daily precipitation variability, where the averaged value from each dataset over grid cells collocated with ground-based observations is marked by a distinct symbol. Shading colour in (a) defines the colour scale used in the spatial distributions of harmonic balance from (b) Multi-Satellite and (c) Multi-Model mean.

To characterize the temporal frequency distribution of regular sub-daily precipitation variability (${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$), the relative contributions of various harmonic components are examined via power spectrum analysis. Ground-based observations show a 1st harmonic (24 h period) fraction of 79.7 %, a 2nd harmonic fraction of 15.1 %, and higher-order contributions of 5.0 % (Fig. 8a). Satellite-based products generally overestimate the contribution of the 1st harmonic relative to observations (MS: 84.1 %), with CMORPH (86.0 %) and TMPA (82.9 %) exhibiting a pronounced overestimation alongside a corresponding underestimation of the 2nd harmonic. Among individual satellite products, MSWEP (78.2 %) captures the observed 1st harmonic fraction well, while also demonstrating superior performance in representing mean climatology (Dirmeyer et al., 2018), suggesting that its multi-source merging approach provides reliable performance across multiple timescales. Within reanalysis products, ERA5 exhibits the lowest 1st harmonic dominance (68.7 %) and the highest 2nd harmonic fraction (25.2 %) of all evaluated datasets, whereas JRA-3Q (87.9 %) presents the most substantial overestimation of the 1st harmonic. While the MM mean (79.4 %) appears numerically closest to observations, this results from the mutual cancellation of ERA5's 2nd harmonic excess and JRA-3Q's 1st harmonic overestimation; therefore, it reflects a statistical artifact rather than a physically consistent reproduction of the observed harmonic balance. The spatial difference against ground-based observations supports this interpretation. Over North America both product suites exceed the observations over the east, but they differ inland, where the MS remains close to or slightly above the observations across the central plains while the MM falls below them. Over Europe the MS differences are small and of mixed sign whereas the MM exceeds the observations over most of the domain, and over East Asia both suites exceed the observations, more so in the MM (Fig. 6d, e). The near-observational domain mean of the MM in Fig. 8a therefore arises from regional contributions of opposite sign rather than from agreement at individual locations.

The spatial distribution of harmonic amplitude ratios in MS (Fig. 8b) indicates that the 1st harmonic dominates over major continental interiors, including the North American Great Plains, India, central Africa, and southeastern China, as well as over midlatitude oceans and continental eastern coasts. Over high-latitude oceans and some subtropical arid regions, no harmonic component exhibits clear dominance, producing irregular patterns that likely reflect insufficient precipitation frequency for stable harmonic decomposition. In contrast, reanalysis products show smaller 1st harmonic amplitudes over major continental interiors compared to the satellite-based retrievals, with correspondingly larger contributions from the 2nd harmonics (Fig. 8c), reflecting the inability of reanalysis convective parameterization to reproduce the observed thermally driven diurnal cycle (Dai, 2024; Tang et al., 2021). The amplified semi-diurnal component in reanalysis products does not reflect a superior representation of physical processes; instead, it indicates the spurious generation of semi-diurnal precipitation variability inherent in model parameterization schemes. Over most ocean basins, the MM exhibits broadly uniform 1st harmonic dominance, including tropical and subtropical regions where the MS displays heterogeneous harmonic compositions. At high latitudes, the MM maintains a comparatively uniform 1st harmonic dominance, contrasting sharply with the spatially incoherent harmonic patterns seen in the MS. This discrepancy between MS and MM is particularly evident over the western Pacific Monsoon Trough, where the 2nd harmonic contribution in the MM exceeds 30 %, whereas the MS exhibits spatially noisy patterns that preclude a clear assessment of harmonic dominance in this region.

Figure 9 Global maps of diurnal phase (colour) and amplitude (colour intensity) of precipitation from (a) MS and (b) MM. Backward slashed area denotes the low precipitation amplitude mask, applied where the daily precipitation amount is less than 0.275 mm or normalized amplitude of diurnal cycle < 30 % (Watters et al., 2021), and forward slash denotes the harmonic dominance mask, applied where rat₁ is less than 50 %. (c, d) Differences in diurnal phase from the ground-based observations for (c) MS and (d) MM, with the mean difference and the spatial RMSE given in each panel. The corresponding maps for the individual datasets are shown in Fig. S3.

4.2 Diurnal phase and amplitude across datasets

The MS captures a pronounced land–ocean contrast in the diurnal phase of precipitation (Fig. 9a). Across most tropical, subtropical, and mid-latitude land regions, precipitation predominantly peaks during the late afternoon to evening (15:00–21:00 LST). In contrast, nocturnal to early-morning maxima (21:00–06:00 LST) are evident over the US GP, parts of the Amazon Basin, central Africa, and the eastern periphery of the Tibetan Plateau. These patterns are underpinned by nocturnal precipitation regimes associated with LLJ activity and propagating MCSs. Over the oceans, precipitation manifests a distinctly different temporal signature, with peaks occurring during the early morning (03:00–09:00 LST), consistent with the nocturnal destabilization of the marine boundary layer driven by radiative cooling aloft (Gray and Jacobson, 1977). Furthermore, the diurnal amplitude exhibits considerable spatial heterogeneity, with the most robust variations observed over land, while oceanic regions are characterized by damped amplitudes due to minimal surface diurnal temperature fluctuations. Along coastal regions, intermediate phases (06:00–12:00 LST) emerge between the land and oceanic regimes, reflecting the offshore propagation of convective systems initiated over land, which is particularly pronounced over the Maritime Continent, western India, West Africa, and the Atlantic coast of South America. The individual satellite products share this structure but differ in detail (Figs. S3a–e). CMORPH and IMERG are the most similar because they share a morphing framework, TMPA carries larger amplitudes, and MSWEP is smoother and peaks near 15:00 LST, inheriting the early timing of ERA5.

The MM peaks earlier than the MS over land while remaining close to it over the ocean, where early-morning peaks at 03:00–06:00 LST align well with the observations (Fig. 9b). It has been commonly reported in evaluations of CMIP models and reanalysis datasets (Tang et al., 2021; Watters et al., 2021), and is associated with convective parameterization schemes triggering precipitation too rapidly in response to solar heating over land, whereas radiatively driven oceanic precipitation timing is more realistically reproduced (Dai and Trenberth, 2004). The nocturnal maxima that the MS resolve over the US GP and the eastern periphery of the Tibetan Plateau are weaker and less spatially coherent in the MM, and the diurnal amplitude over land is lower than in the satellite-based datasets throughout. The phase lead over land varies among the reanalyses (Fig. S3f–h). JRA-3Q remains closest to the satellite products peaking at 14:00–18:00 LST over most continental regions and partially reproducing the nocturnal propagation over the US GP and the North China Plain. MERRA-2 peaks uniformly at 09:00–12:00 LST and resolves neither the nocturnal maxima over mountains nor the afternoon peaks over continental interiors.

Comparing the two suites against the ground-based observations places these phase structures on a common reference (Fig. 9c, d). MS departs from the observations less than MM across the validation domain, with spatial RMSE values of 3.15 h against 4.42 h, and the differences are predominantly positive in MS and negative in MM. The corresponding differences for the individual products are given in Fig. S4. Over the CONUS the MS differences remain within about 2 h across most of the domain and are smallest over the US GP, whereas MM is earlier than the observations by 6 h or more over the upper Midwest. A band of large positive differences in MM runs alongside this area from the southern US GP toward the northeast. MM places its peak in the afternoon across this corridor while the observations peak after midnight, so the two are separated by close to half a day and the wrapped difference changes sign between neighbouring cells. The band therefore marks where MM misses the nocturnal propagation rather than where it peaks late. Over Europe the MS differences stay within about 2 h while MM is earlier by 4 to 6 h over Germany and Poland. East Asia carries the largest MS differences, reaching 6 h or more over northwestern China, where the cells concerned are those hatched in Fig. 9a and b for weak diurnal variability and low precipitation amounts. The harmonic balance over the same cells places a large share of ${\mathit{\sigma }}_{\mathrm{reg}}^{\mathrm{2}}$ in the second and higher-order components (Fig. 8b, c), so a single diurnal harmonic does not describe the sub-daily cycle there and the phase retrieved from it is poorly constrained in both suites. Outside the hatched area, MS is later than the observations by 2 to 4 h over eastern China, so its positive differences concentrate in East Asia while those of MM are negative over almost the whole domain.

Figure 10 Spatially averaged diurnal phase (azimuth) and amplitude (radius) of the precipitation from ground-based observation (black), satellite (blue), and reanalysis (red) datasets over (a) the globe, (b) Contiguous United States (CONUS), (c) Europe (EU), and (d) East Asia (EA), for grid points where the observed rat₁ ≥ 50 %. The spatial correlation of diurnal phase (green solid) and amplitude (red dashed) from gridded products with the observation is depicted in (e–h). The corresponding statistics without this criterion are shown in Fig. S6.

The quantitative analysis of diurnal phase and amplitude corroborates the inter-dataset discrepancies identified qualitatively in Fig. 9. While ground-based observations indicate a peak phase at approximately 17:00 LST, most satellite products, despite capturing the general timing, exhibit a systematic temporal lag of approximately 1–2 h (Fig. 10a). This systematic delay is consistent with the indirect detection characteristics inherent in satellite precipitation retrieval algorithms. PMW sensors detect hydrometeors at ice-scattering altitudes in the upper troposphere rather than at the surface, and preferentially capture CP more effectively than stratiform precipitation (Ebert et al., 2007; Hayden and Liu, 2021). Furthermore, the extensive anvil regions of organized convective systems can continue expanding after surface precipitation reaches its maximum and persist even after the convective core dissipates, such that PMW sensors still detect strong ice-scattering signals during phases when surface precipitation has already weakened. The same indirect detection affects IR-based retrievals, which infer precipitation from cloud-top brightness temperatures that are only weakly related to surface rainfall and tend to place the diurnal maximum later; this is particularly relevant for TMPA, which fills the gaps between PMW overpasses with calibrated IR (Huffman et al., 2007). These detection characteristics cannot be separated from the effects of microwave sampling and temporal interpolation in the present analysis.

The temporal interpolation alone does not account for the delay: TMPA exhibits comparable phasing to the other satellite products despite filling the intervals between overpasses with calibrated infrared estimates rather than propagating retrievals in time, and a comparable lag also appears in CMORPH, whose forward- and backward-propagated retrievals are weighted by elapsed propagation time without preference for either direction (Xie et al., 2017). Backward propagation brings the phase closer to the observations, since the IMERG Early (IMERG-E) Run, which propagates forward only, peaks later than the Late (IMERG-L) and Final (IMERG-F) Runs, which propagate in both directions, although a positive offset remains in all three (Fig. S5). The contribution of sampling density cannot be constrained in the same way, since IMERG, CMORPH, and GSMaP draw their retrievals from the same satellite constellation, and a sampling effect would act on all three products alike. Two products depart from this common behaviour. TMPA overestimates the amplitude, which is a direct consequence of the elevated sub-daily variance identified in Fig. 4a. MSWEP displays a phase of approximately 16:00 LST, roughly 1 h earlier than observations. Reanalysis products exhibit premature phase biases and attenuated amplitudes, with peak precipitation timing at 16:00 LST in JRA-3Q, 14:00 LST in ERA5, and 12:00 LST in MERRA-2.

The diurnal cycle of precipitation exhibits distinct regional characteristics. Over the CONUS, where the peak diurnal phase is observed at 18:00 LST, satellite products generally capture the observed timing but tend to underestimate the robust diurnal amplitude (Fig. 10b). One possible contribution to this deficit lies in the ground-based reference rather than in the products. A grid cell represented by a single station carries the diurnal amplitude of that one record, whereas a cell represented by several carries an average in which maxima occurring at different hours partly cancel, which is closer in character to a gridded estimate, and grid cells over CONUS contain fewer stations than those over the Europe (EU) or the East Asia (EA). The amplitude difference does narrow as the number of contributing stations increases (Fig. S7). But the relationship is weak and a substantial deficit remains where several stations contribute, so the station coverage accounts for only part of the difference.

While most reanalysis products highlight the persistent challenges in simulating diurnal precipitation over this region, NARR demonstrates superior fidelity in reproducing both the observed phase and amplitude. This performance is directly attributable to its 3-hourly assimilation of precipitation observations, thereby underscoring the efficacy of precipitation data assimilation in accurately representing diurnal cycle dynamics.

The observed diurnal phases in EU and EA occur around 16:00–17:00 LST, with magnitudes that are comparable between the two regions yet notably smaller than those over CONUS (Fig. 10c, d). These subdued amplitudes are consistent with the lower rat_subD values observed in the EU (50 %–60 %) and EA (30 %–50 %) compared to CONUS (70 %–80 %). This suggests that synoptic-scale and lower-frequency variability account for a greater proportion of the total precipitation variance in these regions. However, the physical mechanisms underlying the weak diurnal amplitudes differ between the two regions. Over the EU, summertime frontal and orographic precipitation contributes substantially to the total rainfall, whereas the thermally driven convective diurnal signal remains relatively weak. In EA, synoptic-scale systems, such as the Meiyu/Baiu/Changma frontal zones associated with the East Asia summer monsoon, account for a substantial portion of the precipitation variability. Despite the EA exhibiting lower rat_subD than the EU, the two regions manifest similar diurnal amplitudes. This likely reflects the larger ${\mathit{\sigma }}_{\mathrm{total}}^{\mathrm{2}}$ in the EA driven by synoptic environments, while the activation of nocturnal low-level jets partially sustains the diurnal signal even within LSP systems (Yu et al., 2007).

We further evaluate the spatial consistency of the diurnal phase and amplitude between gridded products and observations (Fig. 10e–h). Regarding the diurnal phase, satellite products achieve the strongest correlations over the CONUS (r>0.8), with moderate agreement over EA (r=0.5–0.7) and the weakest performance over the EU. This reflects the distinct spatial structure of the diurnal cycle over CONUS, with an afternoon maximum over the western and eastern parts of the domain and a nocturnal maximum over the central plains, which the satellite retrievals observe directly whereas the reanalyses do not reproduce from surface heating alone. The complex topography and the coexistence of diverse precipitation mechanisms (e.g., frontal, orographic, and convective) in the EU and EA induce higher spatial variability, complicating the representation of the diurnal phase. Reanalysis products attain lower phase correlations than satellite datasets across all four domains. Across the globe, CONUS, and EU, ERA5 shows comparable or higher phase correlations than JRA-3Q, although its domain-averaged peak occurs 2 h earlier, so JRA-3Q remains closer to the observations in the mean phase. An exception occurs in the diurnal amplitude over EA (Fig. 10h), where most datasets (excluding TMPA) yield correlation coefficients approximately 0.2 higher than their respective phase correlations. This indicates that the spatial distribution of amplitude is more consistently reproduced across datasets than the phase.

Figure 11 Diurnal cycle of area-averaged precipitation over land (a–c) and ocean (d–f) grid cells bounded by 60° S–60° N (only grid points after the masking depicted in Fig. 9), for JRA-3Q (a, d), ERA5 (b, e), and MERRA-2 (c, f). Bars show the partitioning of total precipitation (TP) into convective precipitation (CP, orange), and large-scale precipitation (LSP, blue) components, and the black line shows the TP from CMORPH, with the black dot marking its diurnal maximum. Background shaded bars indicate the local time at which each reanalysis reaches its TP maximum.

To identify the physical processes underlying these phase shifts, we examine the respective contributions of the CP and LSP components to the diurnal cycle of TP (Fig. 11). Among the satellite products evaluated above, CMORPH shows the closest agreement with ground-based observations in representing global diurnal characteristics (Fig. 10a), we therefore adopt it as the reference for the grid-scale comparison. Unlike CMORPH, which provides only TP, reanalysis datasets partition precipitation into convective and large-scale components via their respective parameterization schemes. This partitioning enables a detailed investigation into the physical mechanisms responsible for diurnal phase biases. Specifically, this decomposition is instrumental in determining whether these phase errors originate from convective triggering mechanisms, stratiform precipitation processes, or the complex interplay between the two processes.

Over land, CMORPH and JRA-3Q reach their maximum in the 15:00–18:00 LST interval, whereas ERA5 and MERRA-2 peak one interval earlier at 12:00–15:00 LST (Fig. 11a–c). The phase lead of the reanalyses is clearer in the afternoon decline than in the timing of the maximum. Between 18:00 and 21:00 LST, ERA5 and MERRA-2 retain 4.1–4.8 mm d⁻¹ against 8.3 mm d⁻¹ in CMORPH, about half the observed value, whereas JRA-3Q remains comparable at 7.3 mm d⁻¹. Over the ocean, all three reanalyses peak at 03:00–06:00 LST, within the broad early-morning maximum that CMORPH spans from 03:00 to 09:00 LST (Fig. 11d–f), so the shape of the oceanic cycle is reproduced and the phase lead is confined to land. The distributions of the peak hour across grid cells are consistent with these area-averaged results and are given in Fig. S8, where the observed land maximum is broad and carries a secondary late-evening mode over regions influenced by nocturnal MCSs. In every product the diurnal variation of TP over land follows that of CP. Between the daily minimum and maximum of TP, CP varies by 3.1 to 4.4 mm d⁻¹ across the three products, whereas the corresponding change in LSP does not exceed 1.0 mm d⁻¹, so the large-scale component sets the level of the total rather than its phase, and the same alignment between the TP and CP peak phases holds at the grid scale (Fig. S8b–g). ERA5 is the one exception, where CP peaks in the 15:00–18:00 LST but LSP decreases through the afternoon and moves the TP maximum one interval earlier.

The convective share of the total differs among the products. In JRA-3Q, CP accounts for 86 % of TP over both land and ocean, whereas in MERRA-2 it accounts for less than half, so the large-scale component carries most of the precipitation. ERA5 falls between the two and shifts with surface type, from 67 % over land to 56 % over the ocean. The daily-mean LSP in MERRA-2 is 1.8 times that of JRA-3Q over land, which is consistent with the documented tendency of the Relaxed Arakawa-Schubert (RAS) scheme with a Tokioka-type trigger to suppress deep convective activity (Bacmeister and Stephens, 2011; Gelaro et al., 2017) and to shift the burden of precipitation generation toward large-scale condensation (Tao et al., 2024).

The CP-LSP decomposition thus provides direct evidence that premature convective triggering, rather than LSP processes, serves as the primary driver of systematic early peaks in reanalysis TP. This finding corroborates previous studies that identify convective parameterization as the fundamental source of diurnal phase shift (Dirmeyer et al., 2012). To further elucidate the thermodynamic and moisture budget processes underlying these convective phase biases, the following section examines the diurnal evolution of CAPE, MFC, and other related atmospheric variables.

4.3 Physical mechanisms of the diurnal cycle

4.3.1 Moisture-supply regimes over Northern Hemisphere land

To identify the processes behind this phase shift, we examine the diurnal evolution of key variables related to precipitation occurrence, including CAPE, CIN, MFC, and E, on rainy days (daily precipitation ≥ 0.275 mm d⁻¹) over the Northern Hemisphere land (0°–60° N) during the boreal summer (Fig. 12). The variables are averaged over the grid cells of each regime, so they are used for comparison among products rather than as absolute measures of instability. The separation follows the climatological daily-mean ratio MFC $/$ (MFC +E), averaged across the three reanalyses, at a threshold of 0.6 (Fig. S9). Cells above the threshold draw their moisture mainly from horizontal convergence and cells below it from local surface evaporation, so the two regimes correspond to precipitation that is sustained by large-scale transport and to precipitation that is sustained by the surface energy balance.

Figure 12 Diurnal cycle of area-averaged precipitation and its relevant variables for rainy days (daily precipitation ≥ 0.275 mm d⁻¹) over 0°–60° N, separated by the moisture-supply regime into grid cells where the climatological daily-mean ratio MFC $/$ (MFC +E) is (a) at or below 0.6 and (b) above 0.6. Upper panel shows the diurnal cycle of TP (lines) and CP (solid bars) from CMORPH, JRA-3Q, ERA5, and MERRA-2, with circles marking the peak hour of each dataset. Lower panels show diurnal variations of CAPE, CIN, MFC, and E, where each row represents three reanalysis datasets (JRA-3Q, ERA5, and MERRA-2).

The distinction is visible in the observations before any product is considered. In the evaporation-supplied regime, CMORPH reaches its maximum at 15:00–18:00 LST, 3 to 6 h after the surface heating maximum, which is the time required for boundary layer moisture to accumulate, for CIN to erode and for shallow convection to deepen into the free troposphere (Dai et al., 1999; Stirling and Stratton, 2012; Tian and Kuang, 2016). In the convergence-supplied regime, CMORPH reaches its maximum at 21:00–00:00 LST and holds above 7.9 mm d⁻¹ from 15:00 LST through midnight, which reflects the nocturnal precipitation regimes associated with LLJ activity and propagating MCSs. The surface moisture supply does not account for this contrast, since E peaks around 12:00 LST in every product and in both regimes, at 3.9 to 6.6 mm d⁻¹. The two regimes therefore place different demands on the convective schemes, and a single average over the domain conflates them.

In the evaporation-supplied regime, MFC stays close to zero and turns negative during the afternoon in all three products, so the moisture is supplied locally and the phase of the precipitation is governed by the daytime destabilization rather than by the moisture transport. JRA-3Q reaches its maximum at 15:00–18:00 LST, matching the observed phase, whereas ERA5 and MERRA-2 peak 3 h earlier at 12:00–15:00 LST (Fig. 12a). E is alike across the three products, so this difference is not attributable to the surface forcing itself but to how each convective parameterization responds to the resulting destabilization. The evolution of CAPE and CIN describes that response. Their magnitudes differ among the reanalyses despite being computed on the same pressure levels with the same parcel definition and conditional sampling, so the spread reflects differences in the thermodynamic profiles rather than in the diagnostic itself. Consistent with this, the level of free convection (LFC) ordering at 12:00 LST follows the CAPE ordering, with JRA-3Q showing a lower LFC and larger CAPE, indicating that differences in thermodynamic structure contribute to the inter-reanalysis spread in CAPE (Fig. S10).

In JRA-3Q, the two conditions are met at the same hour. CAPE increases gradually after sunrise and reaches its maximum at 15:00 LST (approximately 206 J kg⁻¹) while CIN is weakest at the same hour, and the precipitation maximum follows at 15:00–18:00 LST. CAPE reaches three quarters of that maximum by 09:00 LST, yet the precipitation does not respond while the residual CIN continues to constrain convective development. This behaviour is consistent with the dynamic CAPE generation rate (DCAPE)-based trigger, which requires not only that CAPE exceeds a predefined threshold but also that the rate of CAPE generation by large-scale dynamical processes, such as moisture convergence and free-tropospheric lifting, be sufficiently high. By conditioning the onset on both thermodynamic instability and its rate of generation, the trigger holds convection back through the transition in which shallow cumulus, limited by entrainment of dry environmental air, moisten the lower free troposphere and further erode the residual CIN.

In ERA5 the two conditions are separated. CAPE reaches its peak at 12:00 LST, 3 h earlier than in JRA-3Q, while CIN follows a pattern similar to JRA-3Q, weakening from 12:00 LST and reaching its smallest magnitude at 15:00 LST. TP rises with the CAPE maximum and forms a plateau that holds between 12:00 and 18:00 LST, so precipitation is generated before CIN is fully eroded. CP nonetheless peaks at 15:00–18:00 LST, the same as CMORPH, so the early rise is carried by the total rather than by the convective component. The closure of the modified Tiedtke scheme (Bechtold et al., 2008, 2014) relaxes only the free-tropospheric part of the CAPE, subtracting the contribution generated by boundary layer forcing so that the full instability becomes available to deep convection only at the end of the lower-tropospheric heating and moistening cycle. Where the destabilization originates at the surface, as in this regime, that subtraction holds the convective onset back far enough for CP to peak with the observations, although it does not prevent precipitation from beginning while CIN is still eroding.

MERRA-2 exhibits a different relationship between CAPE and precipitation compared to other products. CAPE reaches its maximum at 18:00 LST, after the precipitation peak at 12:00–15:00 LST, and it continues to increase past the peak while CAPE in JRA-3Q and ERA5 decreases once their precipitation maxima are passed. The convective scheme therefore does not consume the instability at the time precipitation occurs, which is consistent with the RAS convective scheme with Tokioka-type trigger (Bacmeister and Stephens, 2011; Gelaro et al., 2017) to limit the conversion of CAPE into CP, thereby permitting instability to continue building through the afternoon despite the occurrence of precipitation. Instead, the precipitation peak aligns with the maximum phase of E and MFC, both at 12:00 LST, and the strongest midday convergence among the three products, which indicates that MERRA-2 precipitation responds more directly to instantaneous moisture availability than to the cyclic evolution of thermodynamic instability. This behaviour is consistent with the LSP-dominant partitioning (cf., Fig. 11), where suppressed convective activity shifts the burden of precipitation generation toward large-scale condensation processes.

In the convergence-supplied regime, the three products keep the same afternoon timing as in the evaporation-supplied regime, with JRA-3Q peaking at 15:00–18:00 LST and ERA5 and MERRA-2 at 12:00–15:00 LST, while the observed maximum moves to 21:00–00:00 LST (Fig. 12b). None of the schemes respond to the change in the moisture supply, so the phase difference reaches 6 h in JRA-3Q and 9 h in ERA5 and MERRA-2, and the trigger formulations that separate the products in the evaporation-supplied regime no longer distinguish between them. The ingredients for a nocturnal maximum are present. At 21:00 LST, CAPE in all three products remains at three quarters or more of its daily maximum, and MFC in JRA-3Q is largest overnight, more than twice its afternoon value. What changes after the afternoon is CIN, which strengthens in every product between 15:00 LST and midnight. The observed nocturnal maximum requires convection to be initiated against this returning inhibition through the lifting associated with the sustained convergence rather than through surface heating, and none of the three schemes does so. Their precipitation decreases to 4.7–6.4 mm d⁻¹ when CMORPH reaches 8.5 mm d⁻¹, while the afternoon maxima of JRA-3Q and ERA5 exceed the observed daily maximum, so the nocturnal shortage is accompanied by an afternoon surplus.

The two regimes therefore give different answers to the same question. Where the moisture is supplied by local evaporation, the products differ from one another by 3 h and JRA-3Q reproduces the observed phase, so the comparison is informative about the trigger formulations. Where the moisture is supplied by convergence, the three products are indistinguishable from one another and all of them are 6 to 9 h from the observations, because the observed maximum is initiated by sustained lifting rather than by surface heating and none of the schemes represents that pathway. The diurnal phase in the reanalysis products is thus reproduced where precipitation follows the surface energy balance and not where it follows large-scale transport. The same division appears in the phase maps, where MM peaks earlier relative to the observations over the US GP, the region in which the observed maximum occurs at night (Fig. 9d).

4.3.2 Case studies at two sites

The two cases are examined next at the point scale using the ARM SGP site in the US GP, where nocturnal precipitation dominates, and a southeastern US site characterized by afternoon precipitation. The sites and the measurements available at each are summarized in Table 4.

Table 4 Specifications of the two case-study sites examined in Sect. 4.3.2. The observations are available only at the Southern Great Plains site, so CMORPH serves as the precipitation reference at the southeastern United States site. The locations of both sites are marked in Fig. 1.

Figure 13 Same as Fig. 12, but for ARM SGP site, where results from IGRA2 observation and NARR reanalysis are included.

At the SGP site, the diurnal evolution of key variables governing precipitation occurrence is compared with observations (Fig. 13). E shows similar peaks during 09:00–15:00 LST across all products at the SGP site. However, unlike the evaporation-supplied regime where E provides the essential precondition for convective development, nocturnal precipitation dominates at the SGP site and daytime E does not directly contribute to the nocturnal precipitation maximum (cf., Fig. 12). Observations indicate that MFC serves as the primary moisture source for nocturnal precipitation in this region, rather than daytime moisture supply from the land surface. Observed MFC begins to increase after 18:00 LST and reaches peak values up to 20 mm d⁻¹ during 00:00–06:00 LST, coinciding with the precipitation peak. The vertical profiles of moist static energy (MSE) show the same nocturnal phase, with values in the lower troposphere increasing after sunset as the LLJ transports warm and moist air from the Gulf of Mexico (Fig. S11a). The diurnal phase relationship between CAPE and precipitation follows accordingly. In IGRA2 observations, CIN is weakest at 12:00–15:00 LST, as the LFC rises to its highest level while the lifting condensation level (LCL) remains relatively invariant across the diurnal cycle (not shown). Its magnitude varies little through the day, however, so at this site the diurnal phase is governed more by the evolution of CAPE than by the removal of inhibition. Observed CAPE reaches its maximum at 00:00 LST, which is aligned with the nocturnal precipitation peak.

NARR successfully reproduces both the nocturnal CAPE maximum and the MFC peak phase, exhibiting the strongest nocturnal convergence among the products. This results in a precipitation peak at 00:00–03:00 LST that most closely matches the observations. Its MSE bias is confined to 925 hPa and is negligible between 700 and 500 hPa (Fig. S11b), so the vertical moisture distribution associated with the nocturnal LLJ is the most faithfully reproduced of the four products. In contrast, JRA-3Q and ERA5 exhibit weaker nocturnal MFC values compared to NARR. JRA-3Q shows a double-peak structure, with an afternoon peak at 15:00–18:00 LST driven by the CAPE maximum as CIN erodes to its smallest magnitude at 12:00–18:00 LST and a nocturnal peak (03:00–06:00 LST) driven by MFC. However, this nocturnal precipitation phase is delayed by approximately 3 h relative to both the observed peak and its own MFC maximum. ERA5 similarly displays a broad nocturnal peak between 03:00 and 06:00 LST, representing a consistent 3-hour phase delay. Despite a pronounced CAPE maximum near 12:00–15:00 LST, ERA5 does not develop daytime precipitation, instead responding to nocturnal MFC but with a delayed conversion. Both products carry positive MSE biases in the lower troposphere alongside negative biases between 700 and 500 hPa (Figs. S11c, d), so the mid-troposphere moisture transported by the LLJ is underrepresented, and the delayed conversion of the nocturnal convergence into precipitation is consistent with that shortage.

MERRA-2 presents a distinct moisture budget structure. It captures nocturnal MFC but also exhibits an additional strong daytime convergence peak at 12:00 LST and divergence at 18:00 LST, which distinguishes its MFC diurnal evolution from the other products. It exhibits a relatively weak diurnal cycle in CAPE and a mismatch between instability and rainfall. Its MSE bias is positive at every level (Fig. S11e), yet it produces the weakest nocturnal precipitation of the four products, so both the moisture and the instability are present without being converted into precipitation. This is consistent with the lower precipitation efficiency reported for MERRA-2 (Song and Wei, 2021). The persistence of instability without strong convective development indicates that accumulated energy is not efficiently released, consistent with the decoupled CAPE-precipitation relationship identified in the convergence-supplied regime (cf. Fig. 12b).

The three products fail at different points in the same sequence: NARR reproduces both the nocturnal MFC and the vertical moisture structure delivered by the low-level jet, and its precipitation peaks with the observations; JRA-3Q and ERA5 reproduce the nocturnal convergence but convert it into precipitation 3 h late, without the mid-tropospheric moisture that the low-level jet should deliver; and MERRA-2 peaks in the late morning, with moisture present at every level but little of it converted into precipitation at all. Precipitation timing at the SGP site therefore depends on three factors: the representation of nocturnal MFC, the fidelity of mid-tropospheric moisture transport by the low-level jet, and the efficiency with which accumulated instability is converted into CP (Song et al., 2013; Xie and Zhang, 2000). That NARR reproduces all three is consistent with the constraining effect of precipitation assimilation via latent heating nudging. This technique converts observed precipitation into vertical latent heating profiles and incorporates them into the temperature tendency equation, effectively constraining the convective and grid-scale schemes to produce precipitation at observed times (Lin et al., 1999; Mesinger et al., 2006). This result indicates that where convective parameterizations fail to represent the physical pathway connecting moisture transport to precipitation, direct assimilation of precipitation observations can compensate for that structural limitation. Beyond the SGP site, NARR exhibits the smallest phase biases and the highest correlations with observations across the broader CONUS domain (cf. Fig. 10). For representative sites in the southeastern US where IGRA2 observations are unavailable, NARR is therefore used as the reference for diagnosing the diurnal precipitation characteristics of the other reanalysis products.

Figure 14 Same as Fig. 13, but for the southeastern point (33.5° N, 82° W).

The southeastern US point (33.5° N, 82° W), located near Augusta, Georgia, presents a contrasting case to the SGP (Fig. 14). Summertime precipitation at this location is predominantly governed by the influx of abundant lower-tropospheric moisture along the western periphery of the North Atlantic Subtropical High and intense daytime solar heating, which collectively drive CP during the afternoon (Diem, 2006; Li et al., 2012). Unlike the SGP site, convection here is initiated by daytime boundary layer thermodynamic instability (Rickenbach et al., 2015, 2020), though orographic effects from the Appalachian Mountains may also contribute as a secondary factor through eastward-propagating convective systems (Parker and Ahijevych, 2007). This site, therefore, provides a complementary framework for evaluating how reanalysis convective parameterizations respond to boundary layer thermodynamic forcing, specifically for diagnosing the roles of E and MFC.

CMORPH reveals a distinct precipitation maximum at 15:00–18:00 LST, characterized by land-surface thermodynamic forcing. Among the reanalysis products, NARR demonstrates the highest fidelity in reproducing both the phase and amplitude of the observed diurnal cycle. While ERA5 captures the correct peak phase, it consistently underestimates the precipitation intensity, indicating a dampened convective response to surface forcing. Conversely, JRA-3Q exhibits a 3 h phase lead with a peak at 12:00–15:00 LST, and MERRA-2 displays a notably weak, plateau-like distribution spanning 09:00 to 15:00 LST.

The timing and intensity of convection are regulated by the diurnal evolution of boundary layer instability and inhibition. In NARR and ERA5, CAPE reaches a nocturnal maximum and an afternoon minimum, establishing an anti-phase relationship with precipitation. This nocturnal accumulation is sustained by the preservation of boundary layer moisture from daytime evaporation and free-tropospheric radiative cooling. Convective initiation occurs only after surface heating erodes CIN to near-zero values between 12:00 and 18:00 LST, facilitating the rapid release of stored energy. In contrast, JRA-3Q utilizes a DCAPE-based triggering mechanism that activates convection prematurely while CAPE is still being generated by surface heating, resulting in a coupled CAPE-precipitation relationship and reduced convective intensity. MERRA-2, characterized by the largest nocturnal CIN magnitude, undergoes an abrupt reduction in inhibition at 09:00 LST, which triggers a premature and temporally dispersed release of instability, thereby explaining its plateau-like precipitation structure.

The moisture budget at this site further highlights the distinction between thermodynamic preconditioning and dynamical feedback. Across all evaluated products, E peaks between 09:00 and 15:00 LST, serving as a critical preconditioning process that moistens the boundary layer and reduces the gap between the LCL and the LFC. While E establishes the necessary thermodynamic environment, MFC reveals the dynamic response to convective development. In NARR, the MFC peak during 12:00–18:00 LST signifies that convection-induced mesoscale convergence acts as a positive feedback to amplify thermodynamically initiated precipitation. This feedback-driven intensification is followed by an abrupt decline in precipitation at 18:00–21:00 LST, contrasting with the more symmetrical development and decay observed in CMORPH. ERA5 exhibits a qualitatively similar structure, but its consistent underestimation of the afternoon precipitation amplitude is primarily tied to a weaker MFC response. Furthermore, the rapid decay of precipitation in ERA5 after 18:00 LST is driven by its relatively large CIN magnitude, which reflects a faster post-convective stabilization of the lower troposphere compared to other products. On the other hand, JRA-3Q displays prolonged MFC even after the cessation of E during 18:00–00:00 LST. This sustained moisture convergence, combined with its notably weak CIN, facilitates a gradual rather than sharp decay in precipitation throughout the late evening hours.

5 Conclusions

This study evaluates sub-daily precipitation variability and diurnal cycle characteristics across five satellite-based products (IMERG, TMPA, CMORPH, GSMaP, and MSWEP) and four reanalysis datasets (ERA5, JRA-3Q, MERRA-2, and NARR), benchmarked against 6824 ground-based observations over the Northern Hemisphere during boreal summer. By quantifying sub-daily variance fractions, decomposing harmonic structures, and diagnosing the physical mechanisms underlying diurnal phase shifts through atmospheric moisture budget analysis complemented by radiosonde observations, this study advances beyond previous intercomparisons confined to satellite products or limited observation sites.

Ground-based observations establish that sub-daily fluctuations account for 65.6 % of total precipitation variance, and that the irregular component associated with convective intermittency and episodic MCS occurrence dominates (approximately 98 %) over the repeatable diurnal signal. All evaluated datasets underestimate this fraction, most of all the reanalysis products, whose parameterization schemes suppress intense sub-daily events while producing light rainfall too often. Within the regular component, the 1st harmonic dominates across all products, but the balance among harmonics is not reproduced, since ERA5 places too much variance in the semi-diurnal component while JRA-3Q overestimates the 1st harmonic component. The multi-model mean matches the observed first harmonic fraction only because these departures cancel, and a comparable cancellation occurs between regions of opposite sign, so agreement in the domain mean does not indicate agreement at individual locations.

Figure 15 Schematic of the shallow-to-deep transition of convective precipitation over land on a summer afternoon. (a) Surface sensible (SH) and latent (LH) heat fluxes peak near 12:00–13:00 LST as the PBL deepens toward the LCL. (b) Shallow cumulus form once the PBL top exceeds the LCL and moisten the lower free troposphere, limited by the entrainment (ENT) of dry environmental air at cloud top. (c) Congestus clouds deepen further within the moistened layer near 15:00–16:00 LST, with weaker ENT as the clouds broaden. (d) Deep convection develops thereafter, and precipitation, reinforced by downdraft outflow at the surface, reaches its maximum near 17:00–18:00 LST.

The two product classes depart from the observations in opposite directions. Satellite products place the afternoon precipitation peak 1–2 h late, which follows from indirect detection by both PMW sensors, which respond to upper-tropospheric ice-scattering signals rather than surface precipitation, and IR retrievals based on cloud-top temperatures. The delay appears alike in products that fill the gaps between overpasses in different ways, so temporal interpolation does not account for it on its own, although the effect of microwave sampling cannot be separated from the detection physics in this analysis. What separates the satellite products is the design of the correction rather than the sensors. CMORPH agrees most closely with the observations in both diurnal phase and amplitude, which reflects a framework that preserves timing through morphing-based propagation and bias correction approach with PDF matching that leaves hourly magnitudes unconstrained. GSMaP shows a reversed land-ocean contrast in sub-daily variance, plausibly associated with its intensity-dependent retrieval bias over land, namely the overestimation of light precipitation frequency and the underestimation of heavy precipitation intensity, then with its daily-scale gauge adjustment. MSWEP inherits an early phase shift from its use of ERA5 reanalysis, illustrating that the choice of reanalysis input can fundamentally compromise the diurnal fidelity of merged satellite products.

Reanalysis products exhibit a phase 3–6 h early over land, from JRA-3Q (16:00 LST) to ERA5 (14:00 LST) and MERRA-2 (12:00 LST). They also underestimate the absolute magnitude of ${\mathit{\sigma }}_{\mathrm{1}\mathrm{st}}^{\mathrm{2}}$, so the deficiencies extend beyond timing to the intensity. E peaks at the same hour in all three datasets, so the difference lies in how each convective parameterization responds to the resulting destabilization. Figure 15 shows that difference within the sequence separating the surface flux maximum from the precipitation maximum. After the fluxes peak and the PBL deepens past the LCL (a), shallow cumulus moisten the lower free troposphere for 1 to 2 h, a stage held back by the entrainment of dry environmental air at cloud top (b), and only once that layer is moistened can congestus deepen (c) and precipitation reach its maximum (d). The observed lag is the duration of that intermediate stage, and the reanalysis schemes initiate deep convection without completing it.

JRA-3Q departs least from the observations, consistent with its DCAPE-based trigger, which conditions convective onset on both the CAPE threshold and the rate of CAPE generation. That constraint holds the onset through the intermediate stage, but only where moisture is supplied locally at the surface. Comparing the evaporation- and convergence-supplied precipitation regimes shows that this advantage is regime-dependent: JRA-3Q reproduces the observed phase where precipitation is evaporation-supplied, whereas in the convergence-supplied regime all three reanalyses retain their afternoon maxima while the observations peak at night, and the trigger formulation no longer separates the products. The ARM SGP site is one such convergence-supplied location, where NARR uniquely reproduces the observed nocturnal precipitation maximum through direct assimilation of precipitation observations via latent heating nudging. Improving reanalysis diurnal cycles therefore appears to depend on two things: convective triggers that respond to large-scale forcing and not only to local instability, and wider use of precipitation assimilation where the trigger alone is not enough.

The sub-daily fluctuations examined here make up most of total precipitation variance, and every product evaluated in this study underestimates them, whether the precipitation is retrieved from satellite or generated by model physics. Daily and longer averages, which most precipitation evaluations still rely on, do not resolve this part of the variability, even though it is what matters for land-atmosphere coupling (Seo and Dirmeyer, 2022) and for the convective processes analyzed above. Benchmarking several independently developed satellite and reanalysis products against a common set of ground observations helps identify where each product falls short at this time scale and why.

The findings of this study have implications across the disciplines that rely on accurate sub-daily precipitation information. Errors in diurnal phase and intensity in reanalysis products propagate into the representation of land-atmosphere energy and moisture exchanges, bias hydrological simulations of streamflow, and compromise the evaluation of convective parameterization in climate models. The underestimation of sub-daily variance is not confined to one product class but across both satellite and reanalysis datasets, so the representation of precipitation intermittency remains a gap that neither retrieval design nor data assimilation has closed.

This study is limited to Northern Hemisphere land regions during boreal summer, where long-term sub-daily gauge records are dense enough to support the evaluation. Extending the analysis to the Southern Hemisphere, tropical regions, and additional seasons, where the forcing mechanisms and dominant precipitation regimes differ from those examined here, remains a priority for establishing whether the biases identified here generalize beyond the regime sampled. The evaluation of convection-permitting reanalyses, which explicitly resolve mesoscale convective dynamics rather than relying on parameterization, offers a particularly promising avenue for determining whether higher-fidelity process representation can reduce the phase and intermittency biases identified here. Additionally, the rapid emergence of AI-based forecast models such as ECMWF's Artificial Intelligence Forecasting System (AIFS; Lang et al., 2024) raises the question of whether data-driven approaches inherit the phase biases of their training data or instead learn to reproduce the observed diurnal cycle more faithfully.

Appendix A: List of abbreviations

Abbreviation	Meaning
ARM	Atmospheric Radiation Measurement
ARMBEATM	ARM Best Estimate Atmospheric Measurements
CAPE	convective available potential energy
CIN	convective inhibition
CMFD	China Meteorological Forcing Dataset
CMIP	Coupled Model Intercomparison Project
CMORPH	Climate Prediction Center Morphing Technique
CONUS	Contiguous United States
CP	convective precipitation
CPC	Climate Prediction Center
DCAPE	dynamic CAPE generation rate
EA	East Asia
ECMWF	European Centre for Medium-Range Weather Forecasts
ERA5	ECMWF Reanalysis Version 5
EU	Europe
GP	Great Plains
GPM	Global Precipitation Measurement
GSMaP	Global Satellite Mapping of Precipitation
IGRA2	Integrated Global Radiosonde Archive version 2

IMERG	Integrated Multi-satellite Retrievals for GPM
IR	infrared
JMA	Japan Meteorological Agency
JRA-3Q	Japanese Reanalysis for Three Quarters of a Century
LCL	lifting condensation level
LFC	level of free convection
LLJ	low-level jet
LSP	large-scale precipitation
LST	local solar time
MCS	mesoscale convective system(s)
MERRA-2	Modern-Era Retrospective Analysis for Research and Applications, Version 2
MFC	moisture flux convergence
MM	Multi-Model mean
MS	Multi-Satellite mean
MSE	moist static energy
MSWEP	Multi-Source Weighted-Ensemble Precipitation
NARR	North American Regional Reanalysis
NH	Northern Hemisphere
NOAA	National Oceanic and Atmospheric Administration
PBL	planetary boundary layer
PDF	probability density function
PMW	passive microwave
RMSE	root-mean-square error
SGP	Southern Great Plains
TMPA	TRMM Multi-satellite Precipitation Analysis
TP	total precipitation
TRMM	Tropical Rainfall Measuring Mission
UTC	Coordinated Universal Time
VARANAL	ARM Variational Analysis