A Reanalysis Precipitation Integration Method Utilizing the Generalized Three-Cornered Hat Approach and High-Resolution, Gauge-Based Datasets

Abstract

The development of high-precision, long-term, hourly-scale precipitation data is essential for understanding extreme precipitation events. Reanalysis systems are particularly promising for this type of research due to their long-term observations and wide spatial coverage. This study aims to construct a more robust precipitation dataset by integrating three widely-used reanalysis precipitation estimates: Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA2), Climate Forecast System Reanalysis (CFSR), and European Centre for Medium-Range Weather Forecasts Reanalysis v5 (ERA5). A novel integration method based on the generalized three-cornered hat (TCH) approach is employed to quantify uncertainties in these products. To enhance accuracy, the high-density daily precipitation data from the Asian Precipitation-Highly-Resolved Observation Data Integration Towards Evaluation (APHRODITE) dataset is used for correction. Results show that the TCH method effectively identifies seasonal and spatial uncertainties across the products. The TCH-weighted product (TW), calculated using signal-to-noise ratio weighting, outperforms the original reanalysis datasets across various watersheds and seasons. After correction with APHRODITE data, the enhanced integrated product (ATW) significantly improves accuracy, making it more suitable for extreme precipitation event analysis. Quantile mapping was applied to assess the ability of TW and ATW to represent extreme precipitation. Both products showed improved accuracy in regional average precipitation, with ATW demonstrating superior improvement. This integration method provides a robust approach for refining reanalysis precipitation datasets, contributing to more reliable hydrological and climate studies.

1. Introduction

Precipitation is a crucial factor in the water cycle process, and improving the spatiotemporal distribution of precipitation information is essential for enhancing the accuracy of hydrological simulations [1,2,3]. In the context of climate change, the increasing frequency of extreme weather events has increased the demand for high-quality precipitation data [4,5]. Analyzing trends in extreme meteorological and hydrological events, as well as establishing predictive models, necessitate more accurate and comprehensive precipitation datasets to support modeling and statistical analysis [6,7].

Although ground-based rain gauge measurements are the most accurate among precipitation observations, they are spatially and temporally discontinuous, and the gauge stations require regular manual maintenance to ensure data quality. This limitation in station coverage, especially in remote or mountainous regions, further complicates the ability to capture detailed precipitation patterns over large areas [8,9,10]. Currently, with the continuous advancements in satellite remote sensing, atmospheric data assimilation techniques, and data availability, gridded precipitation products (GPPs) provide extensive spatial coverage, high temporal resolution, and real-time accessibility, and are widely available as open-source resources [10,11,12]. These products have greatly supplemented ground-based meteorological stations.

However, GPPs are not directly obtained from observations. For example, satellite precipitation estimations (SPEs) are derived from the remote sensing of rain clouds, and the process of translating these data into ground-level precipitation introduces considerable uncertainties [11,13,14,15]. Similarly, reanalysis precipitation estimations (RPEs) are primarily driven by large-scale meteorological numerical models, and the downscaling process is influenced by the parameterization of various factors [11,16,17]. The accurate representation of convective weather events continues to be a challenge in reanalysis systems [18,19]. To enhance the applicability of these precipitation products, it is crucial to identify their errors and uncertainties and to understand the sources of these errors. This understanding is essential for improving the reliability and accuracy of their usage.

In SPEs, uncertainties often arise from the input data used. High-latitude and high-altitude regions have consistently posed challenges for accurate observation [20]. Studies have shown that as the frequency of precipitation events increases, the ability of remote sensing to capture extreme precipitation diminishes. For events detected by remote sensing satellites, SPEs tend to overestimate light rain while underestimating moderate and heavy rainfall [21,22]. The global precipitation measurement (GPM), as the successor to the Tropical Rainfall Measuring Mission (TRMM), has been credited with improved accuracy in detecting light rainfall through the upgraded advanced microwave (AMW) imager [23]. Apart from improving the accuracy of observational instruments, remote sensing algorithms play a significant role in the uncertainty of results. For instance, while VIS/IR sensors can capture cloud-top temperatures both day and night, and many algorithms estimate precipitation particles based on these data, this indirect observation introduces substantial uncertainty [24]. The PERSIANN system, which relies primarily on infrared data for rainfall estimation, is often reported to be less accurate compared to other satellite products [10].

For RPEs, the observational data inputs for reanalysis systems are irregular in both space and time. The data assimilation process combines available observations with numerical weather prediction (NWP) forecasts to generate a consistent gridded dataset. As a result, while the output data are guided by observations, the physical processes within meteorological models and inherent uncertainties still contribute to the final product’s uncertainty [17]. Betts et al. summarized the strengths, weaknesses, and hydrological impacts of reanalysis precipitation estimates (RPEs) [25]. Precipitation is a key component of both the water and energy cycles but is also strongly tied to simulated physical parameters. Newman et al. found internal consistency among three different reanalyses in precipitation, outgoing longwave radiation, and upper-level divergence, but the agreement between the reanalyses was generally low [26]. Chen et al. (2008) showed that the Hadley circulation in the European Centre for Medium-Range Weather Forecasts (ECMWF’s) 40-year reanalysis (ERA5-40) had significant temporal shifts, but these shifts were likely more related to false trends in precipitation (and the corresponding latent heat) rather than actual atmospheric changes [27]. This suggests that precipitation, as a crucial element of both the energy and water cycles as well as the dynamic circulation, should be a key indicator of the quality of reanalysis for climate studies. Reanalysis precipitation is closely linked not only to the system’s physical processes but also to data assimilation. Although most current global reanalyses do not directly assimilate precipitation—treating it instead as an output variable (except for datasets using the latest 4D-VAR technology)—the data assimilation process still significantly influences the reanalysis precipitation output [17].

There are distinct spatiotemporal differences between SPE and RPE products. RPEs generally exhibit lower accuracy than SPEs over annual scales, particularly in their ability to reflect both rainy and non-rainy events; coarse resolution also hampers local climate representation. However, RPEs often provide decades of model output, which is crucial for analyzing long-term climate change, especially when high-resolution observational data are scarce. To address uncertainties stemming from indirect observations or data assimilation, major research institutions currently focus on correcting reanalysis data using ground-based observational datasets [17,23]. This approach has proven to be highly effective, although it is limited by the delay in releasing ground-truth datasets, meaning that corrected products typically lag by a few days to several months. High-density regional rain gauge networks are also considered highly valuable. Additionally, to correct errors in extreme value estimations, some studies have employed distribution-based mapping methods, yielding promising results [28]. To improve observation and forecasting of the hydrometeorological environment in Asia, the APHRODITE (Asian Precipitation-Highly-Resolved Observation Data Integration Towards Evaluation) project aims to develop an advanced gridded precipitation dataset with a spatial resolution of 0.25° × 0.25° (Lon × Lat) and daily temporal resolution, based on extensive ground observation data across the region [29]. Since its release, the daily precipitation dataset from APHRODITE has been widely applied, particularly in mountainous areas like the Tibetan Plateau. Numerous studies have focused on integrating multiple rainfall products to create improved precipitation datasets [30].

In summary, GPPs and ground-based observations each offer unique advantages and limitations, making them complementary for assessing precipitation [3,31]. Due to the limited temporal scope of satellite remote sensing and the insufficient long-term coverage of ground-based precipitation measurements, reanalysis products serve as vital resources for studies on global climate change. Consequently, this study focuses on improving RPE products by proposing an integration method to obtain more reliable, long-term precipitation data suitable for hydrological and climate research.

2. Materials and Methods

2.1. Study Area

China, with its vast territory and diverse geographical features, exhibits varying topography, including plateaus, plains, and flatlands, with an elevation gradient sloping from west to east. Geographically, China can be divided into three topographic steps. The country experiences a range of climates, which are categorized into three natural zones: the eastern monsoon region, the arid and semi-arid northwest region, and the high-altitude mountain climate region.

To investigate the accuracy of different precipitation products across various regions, this paper divides China into nine river basins based on its hydrological characteristics. The nine basins, from north to south, are the Song-Liao River Basin, the Inland River Basin, the Hai River Basin, the Yellow River Basin, the Huai River Basin, the Yangtze River Basin, the Southwest River Basin, the Southeast River Basin, and the Pearl River Basin.

Figure 1a illustrates the elevation and the annual average precipitation across these nine basins, while Figure 1b shows the spatial distribution of the annual mean precipitation from 2008 to 2015. Precipitation in mainland China increases from northwest to southeast, and from west to east. Due to the subtropical monsoon climate, the Pearl River Basin and the Southeast River Basin in the southeastern part of China are the wettest areas. In contrast, the Inland River Basin, which is far from the influence of the monsoon, receives only one-tenth of the average precipitation of the Pearl River and Southwest River basins.

The basin divisions and elevation data were sourced from the Resources and Environmental Sciences and Data Center of the Chinese Academy of Sciences: https://www.resdc.cn/ (accessed on 1 August 2024).

2.2. Data Description

2.2.1. Reanalysis Datasets

Reanalysis products offer long-term observations of multiple climate variables, primarily derived from the output of numerical models and assimilated with multi-source data for correction. The characteristics of reanalysis products make them well-suited for analysis related to climate change, offering promising application prospects. In this paper, three mainstream RPE products are utilized for evaluation and integration.

The selected RPEs are Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA2) from the National Aeronautics and Space Administration (NASA), Climate Forecast System Reanalysis (CFSR) from the United States National Center for Environmental Prediction (NCEP), and European Centre for Medium-Range Weather Forecasts Reanalysis v5 (ERA5) from ECMWF.

MERRA2 is generated based on NASA’s v.5.2.0 Goddard Earth Observing System (GEOS) atmospheric model and data assimilation system [32,33], and provides hourly data starting from 1980, with a spatial resolution of 0.625° × 0.5° (longitude/latitude). CFSR is developed by NCEP and can be categorized into two versions based on the release timeline: CFSv1, which covers data from 1979 to 2010, and CFSv2, which includes data from 2011 to the present. Since both versions are used in the comparative analysis with reference datasets that span both records, they are collectively referred to as CFSR. CFSR is based on a fully coupled ocean–land–atmosphere model and employs numerical weather prediction (NWP) techniques to assimilate and forecast atmospheric states [34,35].

The CFSR reanalysis system first utilizes a high-resolution Global Data Assimilation System (GDAS) to temporally and spatially decompose model precipitation, which has good quality but lower resolution. It then applies a latitude-dependent linear weighting method to integrate these precipitation data, leveraging the advantages of satellite, rain gauge, and model precipitation across different regions, especially in rain gauge-dense areas like the continental United States, Europe, and Australia, giving the maximum weight to rain gauge analyses, thereby significantly improving precipitation accuracy. The next generation, MERRA2, adopts a similar driving approach to enhance rainfall simulation.

ERA5, the latest fifth-generation global atmospheric reanalysis product released by the ECMWF, succeeds ERA-Interim and is currently a leading reanalysis product [36]. A key feature of ERA5 in precipitation data processing is its use of a four-dimensional variational (4D-Var) data assimilation system. Unlike the three-dimensional variational (3D-Var) system employed by MERRA2 and CFSR, where rainfall is treated as a model output variable not directly influenced by observational data, the 4D-Var system incorporates precipitation observations into the assimilation process. The NWPs of MERRA2 or CFSR do not directly assimilate precipitation observations into the atmospheric model. Instead, they use observed precipitation data, such as those from the Climate Prediction Center (CPC) ground gauge dataset or Global Data Assimilation System (GDAS) products or satellite-derived precipitation, to correct the land surface model component. The 4D-Var assimilation system of ERA5 allows the analysis to account for both model predictions and physical processes impacting rainfall, significantly improving the forecasting of precipitation intensity and the positioning of rain bands, which has a notable positive effect on precipitation predictions for mesoscale convective systems.

Basic information on the above RPEs is listed in

Table 1. Summary of basic information of gridded precipitation products used in this study.

Product	Version	Spatial Resolution	Temporal Resolution	Record from	Institute
CMPA	V1.0	0.1° × 0.1°	1 h	2008.01	CMD
APHRODITE	V1.0	0.25° × 0.25°	1 day	1951.01	Yatagai et al. [29]
MERRA2	V2, LAND	0.625° × 0.5°	1 h	1980.01	NASA
CFSR	V1/V2	0.312° × 0.312°	1 h	1980.01	NCEP
ERA5	V5, LAND	0.1° × 0.1°	1 h	1950.01	ECMWF

.

2.2.2. A High-Resolution and Gauge-Based Dataset

As mentioned above, ground gauge-based datasets are commonly used in correction datasets. Although both the CPC and GPCC have developed a series of gauge-based precipitation analysis datasets with continuously improving quality and spatial-temporal resolution, accurately estimating terrestrial precipitation remains challenging due to the limitations of rain gauge networks [30].

To more accurately observe and simulate the hydro-meteorological processes in Asia, the Asian Precipitation-Highly-Resolved Observation Data Integration Towards Evaluation (APHRODITE) project aims to develop a state-of-the-art GPP with a spatial resolution of 0.25° × 0.25° (Longitude × Latitude) and a daily temporal resolution, covering the entire Asian continent. The dataset is based on a large volume of ground-based observational data from across Asia [29]. Since its release, the APHRODITE daily precipitation dataset has been widely applied, particularly in topographically complex areas like the Tibetan Plateau [29]. Numerous studies have focused on integrating multiple rainfall products to obtain improved precipitation datasets [3,30,37,38].

The APHRODITE data are available on the website http://aphrodite.st.hirosaki-u.ac.jp/download/ (accessed on 1 August 2024). Basic information on APRODITE is listed in Table 1.

2.2.3. Benchmark Precipitation Dataset

To compare the performance of various RPEs and the integrated product, the China Hourly Merged Precipitation Analysis (CMPA) dataset, which merges data from China’s automatic weather stations (AWS) and CMORPH, was used as the benchmark precipitation dataset. CMPA has a spatial resolution of 0.1° × 0.1° (longitude × latitude) and a temporal resolution of 1 h. It combines data from over 30,000 AWSs with CMORPH through a Probability Density Function-Optimal Interpolation (PDF-OI) method, where more than 50% of the CMPA grids contain multiple ground-based automatic stations [39].

As mentioned earlier, the integration process of automatic stations and CMORPH in CMPA may affect the independence of each product’s evaluation. Considering that CMORPH mainly provides spatial interpolation for areas without stations, in order to avoid the influence of the integration process on the ground station dataset, it is generally recommended to use grids with as many stations as possible as the reference precipitation dataset [10,37,40]. According to Tian et al. [41] and Mandapaka et al. [42], spatial mismatches between the grid and the stations within it may introduce uncertainties in the rainfall results. These uncertainties can be minimized when the number of stations within a grid is four or more.

In this study, out of the more than 30,000 AWSs, 2169 are national basic stations. Typically, ground-based automatic stations are built around these basic stations to densify the observation network. Grids containing these basic stations are most likely to have more automatic stations. Therefore, the CMPA grids containing national basic stations were extracted and used as the ground-based observed rainfall data. Due to some basic stations being located within the same grid, a total of 2144 CMPA grids were ultimately selected.

Basic information on CMPA is listed in Table 1. The selected AWSs are shown in Figure 1b.

2.3. Integration Method

2.3.1. Three-Cornered Hat Method

The generalized three-cornered hat (TCH) method is a commonly used approach for assessing the uncertainties of remote sensing products or hydrological models [43,44]. TCH was originally proposed for measuring the instability of clocks. The classical three-cornered hat method evaluates the uncertainty of three observational sequences in the absence of a true value. Later, a study expanded the application scope of TCH by minimizing the global correlation to account for the relationships between datasets, extending it to scenarios involving any number of evaluation objects [45]. There are currently very few studies that have tested its effectiveness in the evaluation of precipitation [38]. In this paper, the uncertainty of each RPE is evaluated using the TCH method, and the weights for integration are determined based on the signal-to-noise ratio. Assuming there are N precipitation products $\left\{P_i,i=1,2,\dots ,N\right\}$ to be evaluated, $P_i$ can be written as follows:

$$P_i=P_{true}+{\epsilon }_i$$

(1)

where $P_{True}$ represents the unknown true precipitation, ${\epsilon }_i$ is the observation error for the ith product sequence. Since the true value is not known in the application of TCH, one of the precipitation product sequences is typically chosen as the reference. The difference between this reference and the other N−1 products is given by the following:

$$D_{iN}=P_i-P_N={\epsilon }_i-{\epsilon }_N,i=1,2,\dots ,N-1\left\{P_i,i=1,2,\dots ,N\right\}$$

(2)

where $P_N$ represents the selected reference precipitation product sequence. According to the derivation by Ferreira [43], the choice of the reference sequence does not affect the calculation results of TCH. Assuming there are M time steps in the current time series, Equation (2) can be expressed as a matrix of size $M\times (N-1)$:

$$D=\left[D_{1N} D_{2N} \cdots D_{(N-1)N}\right]$$

(3)

The covariance of the difference matrix D can be calculated using the following equation:

$$S=\left[{\displaystyle \frac{1}{M-1}}\right]\left[{(D-\overline{D})}^T(D-\overline{D})\right]$$

(4)

where $\overline{D}$ is the mean of matrix D. Here, we introduce an unknown symmetric noise covariance matrix R, the relationship between R and S can be expressed as follows:

$$S=K^T\cdot R\cdot K, K=\left[\begin{array}{l}I\\ -u^T\end{array}\right]$$

(5)

In Equation (5), $u={\left[1\cdots 1\right]}^T$, and I is the identity matrix. A detailed derivation of the above expression can be found in the literature by Premoli and Tavella [45]. Therefore, Equation (5) can be expressed as follows:

$$S=\left[I,-u\right]\left[\widehat{R}r{r}^T{r}_{NN}\right]\left[\begin{array}{l}I\\ -u^T\end{array}\right]$$

(6)

where $\widehat{R}$ is the covariance matrix of size $(N-1)\times (N-1)$, and $r=[r_{1N},r_{2N},\cdots ,r_{N-1,N}]$ represents the group covariance estimate between each sequence and the N-th time series, with $r_{NN}$ denoting the variance of the N-th time series. Once $r$ and $r_{NN}$ are determined, $\widehat{R}$ can be computed using the following equation, and Equation (5) can be solved as follows:

$$\widehat{R}=S-r_{NN}\left[u{u}^T\right]+u{r}^T+r{u}^T$$

(7)

Tavella and Premoli [46] proposed solving for the N-free parameters in $r$ and $r_{NN}$ by imposing the condition of positive definiteness on the covariance matrix (i.e., det(R) > 0).

Premoli and Tavella [45] introduced an optimal criterion to minimize the global correlation between the noise of the time series while maintaining the positive definiteness of R. Building on this, Galindo and Palacio [47], using the Kuhn–Tucker theory, presented an improvement that minimizes the quadratic mean of the covariance of the free parameters as expressed in Equation (8):

$$F(r,r_{NN})={\displaystyle \sum _{i<j}\frac{r_{ij}^2}{L^2}}$$

(8)

and is subject to the following constraints:

$$G(r,r_{NN})=-{\displaystyle \frac{r_{NN}-{\left[r-r_{NN}u\right]}^T\cdot S^{-1}\cdot \left[r-r_{NN}u\right]}{L}}<0$$

(9)

In Equation (8), $L=\sqrt[N-1]{\det (S)}$, following the method described in research by Torcaso et al. [48], the initial conditions for the iterative process are set as follows:

$$\left\{\begin{array}{l}{r}_{iN}^0=0,i<N\\ r_{NN}^0={\left(2\cdot u^T\cdot S^{-1}\cdot u\right)}^{-1}\end{array}\right.$$

(10)

Once the free parameters are determined, the unknown elements in R can be computed using Equation (7).

After estimating the uncertainty of the rescaled data using the TCH method, the relative weights of each product are determined by calculating the signal-to-noise ratio (SNR). A higher SNR indicates that the product contains more signal and should be assigned a higher weight. In this study, the SNR is calculated by dividing the standard deviation of each product by the uncertainty determined using the TCH method [43].

2.3.2. Integration of Three Reanalysis Precipitation Estimation

This paper aims to obtain a more accurate precipitation dataset based on MERRA2, CFSR, ERA5, and APHRODITE. The improved approach and steps are as follows:

(1) Resample MERRA2, CFSR, ERA5, and APHRODITE to a consistent spatial resolution of 0.1° × 0.1° (Longitude × Latitude). Additionally, CMPA is a product that integrates CMORPH and high-density AWS gauges, as shown in Figure 1b, and many areas in the northwestern mountainous regions of China lack AWS coverage. During the evaluation phase, to avoid the influence of the fused products on spatial rainfall, this study only selects grids with AWS coverage from the resampled products for analysis.

(2) To integrate the three RPE products, calculate the uncertainty and SNR for each product based on the TCH. The weights for the three products will be determined according to the SNR, and a weighted average of the precipitation for each grid and time period will yield the TCH-weighted product, abbreviated as TW.

(3) Correct the results of TW based on APHRODITE:

(3-1) Rain occurrence/non-occurrence correction, where all precipitation for dates with no detected rainfall in APHRODITE is set to zero;

(3-2) Daily rainfall correction, where hourly TW rainfall is aggregated to daily scale (noted as D-TW), and the hourly rainfall of TW is redistributed based on the daily D-TW relative to APHRODITE’s daily rainfall, ensuring the total daily rainfall matches APHRODITE for the corresponding date, resulting in the final product: APHRODITE calibrated TW, referred to as ATW.

The overall computation flow is illustrated in Figure 2. The precipitation products used in this study span from 2008 to 2015. Specifically, rainfall data from 2008 to 2011 (four years) will be used to calculate TCH and weights, while these weights will be applied to the 2012–2015 data to obtain final results for comparison with the CMPA-extracted results.

2.4. Evaluation Metrics

To evaluate each product, this study selects statistical indicators from

Table A1. Statistical metrics.

Evaluation Metric	Value Range	Optimal Value	Formula
RMSE	$[0,+\infty )$	0	$RMSE=\sqrt{{\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{\left(T_i-R_i\right)}^2}}$	(A1)
C-C	$[0,1]$	1	$\mathit{\text{C-C}}={\displaystyle \frac{{\displaystyle \sum _{i=1}^n(R_i-\overline{R})(T_i-\overline{T})}}{\sqrt{{\displaystyle \sum _{i=1}^n{(R-\overline{R})}^2}}\cdot \sqrt{{\displaystyle \sum _{i=1}^n{(T-T)}^2}}}}$	(A2)
PBIAS	$(-\infty ,+\infty )$	0	$PBIAS={\displaystyle \frac{{\displaystyle \sum _{i=1}^n(R_i-T_i)}}{{\displaystyle \sum _{i=1}^n{T}_i}}}\times 100\%$	(A3)
$KGE \prime$	$(-\infty ,1]$	1	$KGE\prime =1-\sqrt{{\left(r-1\right)}^2+{\left(\beta -1\right)}^2+{\left(\gamma -1\right)}^2}$ $\beta ={\displaystyle \frac{{\mu }_{tar}}{{\mu }_{ref}}}$ $\gamma ={\displaystyle \frac{{\mathrm{CV}}_{tar}}{{\mathrm{CV}}_{ref}}}={\displaystyle \frac{{\sigma }_{tar}/{\mu }_{tar}}{{\sigma }_{ref}/{\mu }_{ref}}}$	(A4)
				(A5)
				(A6)
POD	$[0,1]$	1	$POD={\displaystyle \frac{H}{H+M}}$	(A7)
FAR	$[0,1]$	0	$FAR={\displaystyle \frac{F}{H+F}}$	(A8)
CSI	$[0,1]$	1	$CSI={\displaystyle \frac{H}{H+M+F}}$	(A9)

based on existing research for comparative assessment. The root-mean-square error (RMSE) is a common measure of the differences between numerical values, representing the quantitative discrepancy between estimated and observed values, as expressed in Equation (A1) in Table A1. The correlation coefficient (C-C) is one of the most commonly used statistical indicators for assessing precipitation product performance, characterizing the linear relationship between estimates and observations; the Pearson correlation coefficient formula used in this study is given in Equation (A2). RMSE cannot indicate whether the estimates are overestimations or underestimations relative to the measured values; therefore, this study introduces the percentage bias (PBIAS) to describe the direction and magnitude of deviation of estimates relative to observations (Equation (A3)). The Kling–Gupta efficiency (KGE’) is a comprehensive metric that considers the mean, correlation coefficient, and coefficient of variation of estimates [49,50]. In this study, KGE’ is calculated for each evaluated product based on Equations (A4) to (A6), where $\mu$, $\sigma$, and CV represent the mean, variance, and coefficient of variation of the response dataset, respectively. The subscripts “tar” and “ref” denote the target dataset being evaluated and the reference dataset, which in this study corresponds to the observed rainfall obtained from CMPA (the definitions of the subscripts apply throughout this study).

The classification of rain and no-rain events is a critical challenge faced by the inversion algorithms of various precipitation products. In this study, a commonly used binary classification metric method is adopted for evaluation. This method, proposed by Wilks [51], has been widely applied in subsequent research. The primary idea of this approach is to first assume a threshold for rain occurrence. When both the estimated precipitation from a product and the observed precipitation exceeds the threshold, the product is considered to have successfully identified the rain event (Hit alarm, denoted as H). If the estimated precipitation exceeds the threshold but the observed value does not, the product is considered to have issued a false alarm (false alarm, denoted as F). Conversely, if the product’s estimated precipitation is below the threshold while the observed precipitation exceeds it, the product is considered to have missed the rain event (miss alarm, denoted as M). After obtaining the detection results of the precipitation product relative to observed rainfall, its detection accuracy can be evaluated using the probability of detection (POD), false alarm ratio (FAR), and critical success index (CSI), calculated by Equations (A7) to (A9).

2.5. Quantile Mapping

The quantile mapping (QM) correction algorithm is typically used to adjust systematic distribution biases in climate model precipitation outputs [28]. According to Zhang et al. [12], comparing satellite remote sensing and other observed rainfall products, the potential limitation of reanalysis products lies in their poor simulation of extreme values. This study aims to improve the accuracy of extreme precipitation modeling through quantile mapping.

QM is a method that equates the cumulative distribution function (CDF) of observed data to the CDF of simulated data. Using the subscript o to denote observed data, m for simulated data, h for the calibration period (i.e., the period used to obtain the CDF parameters), and p for the mapping period, the QM of simulated data at time t during the mapping period can be expressed by the following equation:

$${\widehat{x}}_{m,p}(t)=F_{o,h}^{-1}\left\{F_{m,h}\left[x_{m,p}(t)\right]\right\}$$

(11)

where $F_{o,h}^{-1}$ is the inverse CDF of the observed data; $F_{m,h}$ is the CDF of the simulated data during the calibration period; $x_{m,p}(t)$ represents the simulated data at timestep t during the mapping period after the QM correction.

3. Results

3.1. Product Uncertainty Evaluation

The uncertainty of MERRA2, CFSR, and ERA5 at a 0.1° × 0.1° resolution for the period from 2008 to 2011 was calculated using the generalized three-cornered hat method and is illustrated in Figure 3. Due to significant differences in accuracy across time and space, this study calculates uncertainty based on seasons and river basins. MERRA2 consistently exhibits lower uncertainty than other products across all four seasons, particularly in the lower reaches of the Yangtze River and in southern China, where uncertainty is highest. MERRA2 reported the lowest rainfall error, and in a previous study by Xu et al. [38] on global monthly rainfall errors, MERRA2 also showed the lowest rainfall error. It is important to note that the magnitude of uncertainty is related to the average rainfall of the target grid, and since true values are unknown, first-order systematic bias cannot be computed; therefore, the uncertainty estimated based on TCH should be less than that estimated from true values.

The spatial distribution of the SNR for each product from 2008 to 2011 is depicted in Figure 4. For each individual product, the SNR is determined as the ratio between the standard deviation of the rainfall time series for each grid and the uncertainty identified in Figure 3. The standard deviation represents the useful information (signal) present in the dataset, while the uncertainty quantifies the level of noise. A higher SNR indicates that the dataset contains a more substantial amount of useful signal compared to noise, thereby suggesting greater accuracy and reliability of the product. Conversely, a lower SNR suggests that noise dominates, reducing confidence in the data’s accuracy and reliability. In cases where the SNR is less than 1, it implies that the noise level surpasses the signal in that grid point, rendering the data less dependable. These regions are marked with hollow circles in Figure 4 to highlight areas where the information-to-noise ratio is less favorable. The results indicate that while MERRA2 generally exhibits the lowest uncertainty, its SNR is lower than that of the other two products across most regions, particularly during winter. However, in the northern region, and specifically in the Song-Liao River Basin, MERRA2 demonstrates a higher SNR compared to the other products during the spring and autumn seasons. In the southwestern river basins and the upper reaches of the Yangtze River during winter, MERRA2 exhibits a large area with an SNR of less than 1. In contrast, ERA5 shows the highest SNR, particularly in the Song-Liao River Basin during winter, where high signal values are evident.

As previously analyzed, the low uncertainty of MERRA2 may be related to its coarse resolution, which allows it to provide more reliable estimates over large areas. However, this coarse resolution limits its ability to capture local climatic characteristics and small-scale rainfall events. Therefore, the following sections aim to combine coarse-resolution products with higher-resolution products to obtain more accurate rainfall estimates.

Based on the SNR results in Figure 4, this paper calculates integration weights for each RPE product, with each grid point having its own set of weights. To illustrate the seasonal variations in these weights, Figure 5 presents the seasonal weighting of individual RPE rainfall products derived from the SNR values for the period 2008 to 2011. These weights are then used to perform a weighted average of MERRA2, CFSR, and ERA5, resulting in the TW product. The weights reflect the relative uncertainties of each rainfall product, similar to the information shown in Figure 4. Across mainland China, MERRA2 has the highest weight in summer and autumn, while ERA5 carries a higher weight in spring and winter. As shown in Figure 4, MERRA2 demonstrates a higher SNR in the mid-to-low latitude regions during summer, while the higher winter weighting of ERA5 is attributed to its relatively high SNR in the Song-Liao River Basin.

3.2. Daily Rainfall Correction

Gridded rainfall products derived from high-density ground stations have played a crucial role in correcting other rainfall products [30,52]. Based on the weights obtained in the previous section, the combined weighted (TW) rainfall sequence was derived. Assuming the rainfall sequence from APHRODITE, a gridded product based on high-density ground station data, represents the true values, the daily rainfall sequences from MERRA2, CFSR, and ERA5 were aggregated to the daily scale. The ratios between these rainfall sequences and APHRODITE daily rainfall were calculated, as shown in Figure 6.

The results show that the correction factors for RPEs exhibit significant seasonality, with daily average rainfall generally higher than APHRODITE observations. In mid- and high-latitude regions, the correction factor for MERRA2 remains consistently stable across different basins and seasons, whereas in mid- and low-latitude regions, including the Yangtze River Basin, southwestern rivers, southeastern rivers, and the Pearl River Basin, all original RPE products demonstrate instability.

3.3. General Metrics

Based on the results of TW and Figure 6, the hourly rainfall product ATW, obtained through TCH integration and APHRODITE correction for the period 2012–2015, is generated. This section compares the accuracy of MERRA2, CFSR, ERA5, TW, ATW, and observed rainfall extracted from CMPA using the statistical evaluation metrics listed in Table A1, to assess the improvement achieved. Figure 7 shows the monthly average precipitation of three RPEs TW and ATW in mainland China.

The monthly average results indicate that ATW values are closer to the observed monthly rainfall sequence. However, according to previous studies, when the data were aggregated on a monthly scale, the differences between various reanalysis products were minimal [10,12,20]. Therefore, Figure 8, Figure 9 and Figure 10 display seasonal evaluation metrics such as C-C, $KGE’s$, and CSI. The spatial distribution of KGE in each season is shown in Figure 11. TW shows improved accuracy compared to original RPE products, but ATW significantly outperforms the others, especially in the improvement of C-C.

As shown in Figure 8, Figure 9 and Figure 10, TW enhances the average metrics, but its optimization for extremes and the issue of rain/no-rain events is uncertain. In Figure 10, TW’s CSI does not show a clear improvement over other RPE products.

To further analyze the improvements in ATW across different regions, accuracy metrics for mainland China, Songhua River Basin, Inland River Basin, Southwest Rivers Basin, Yangtze River Basin, and Pearl River Basin are shown seasonally in Figures S1–S4. The scatter plots demonstrate that TW and ATW data points are closely clustered around the y = x line, indicating high accuracy for both. The calculated accuracy metrics reveal that ATW provides a noticeable improvement over TW.

As noted earlier, due to insufficient representation in sub-daily precipitation, RPE products score significantly lower $KGE\prime$ on the hourly scale, particularly during the winter season, with occasional negative $KGE\prime$ values. TW did not show significant improvements in negative $KGE\prime$, but ATW displayed substantial enhancements. Comparing the average $KGE\prime$ values for mainland China, ATW outperformed TW by approximately 0.15 across all seasons.

3.4. Extreme Event Correction

This section evaluates the interpolation between ATW and observed rainfall sequences based on the quantile mapping (QM) approach and explores the potential for improving ATW through post-processing methods.

Firstly, the average rainfall for each time period was calculated across stations in various regions, and the percentiles were plotted as shown in Figure 12. The results indicate that ATW tends to underestimate rainfall above the 50th percentile across different regions.

Secondly, for selected regions (e.g., mainland China, Song-Liao River Basin, Inland River Basin, Yangtze River Basin, Southwest River Basin, and Pearl River Basin, which either exhibit significant dry spells or are difficult to predict), the hourly rainfall from all stations was combined into a single sequence. The percentile process was then computed to represent ATW’s ability to capture extreme rainfall at individual stations within these regions, as depicted in Figure 1. The results show that ATW tends to underestimate rainfall above the 30th percentile.

To address this, a quantile mapping (QM) approach was applied to correct the values at each percentile. As shown in Figure 12 and Figure 13, the QM-adjusted ATW (QM-ATW) provides improved accuracy in terms of regional average rainfall. However, in regions with higher rainfall, such as the Yangtze River Basin and Pearl River Basin, the improvement is only marginal. While QM offers a useful approach for adjusting extreme values in RPE and related rainfall products, the maximum rainfall after adjustment may still exceed the corrected value.

4. Discussion

Obtaining long-term, stable, high-resolution precipitation data is crucial for climate change analysis and hydrological model calibration. While long-term reanalysis datasets present a promising solution, evaluations have identified several limitations in reanalysis precipitation estimates (RPEs) [17,53]. Integrating different RPE products can help address model-specific shortcomings and has wide applications in hydrology and meteorology [37,38,54]. This study utilizes the TCH method, alongside the high-resolution ground gauge-based precipitation dataset APHRODITE, to enhance the accuracy of RPE products. By applying these methodologies, the study further assesses the improved product’s capacity to represent extreme precipitation events.

This study integrates three reanalysis precipitation products—MERRA2, CFAR, and ERA5. The TCH method was first applied to estimate the uncertainties of the three RPEs, which were then used to calculate SNR-based weights. The resulting TCH-integrated product, referred to as TW, generally demonstrates a higher correlation and lower RMSE than the original products. However, as noted in prior research, TCH may underestimate uncertainty due to representation errors [55,56] and is limited in identifying specific causes of uncertainty, necessitating further statistical analysis. Additionally, simply integrating reanalysis products via TCH cannot fully address the inherent limitations of RPEs. To enhance TW, the high-resolution APHRODITE ground gauge-based precipitation dataset was incorporated for correction. Results indicate that including high-quality observational data significantly improves product accuracy, with APHRODITE effectively addressing the issue of rainfall occurrence/non-occurrence. General evaluation metrics further show that ATW achieves more consistent improvements than TW, markedly reducing RPE product errors across all basins and seasons.

Extreme precipitation events can lead to natural disasters such as floods, urban waterlogging, and flash floods [57,58], and are also linked to the health of terrestrial ecosystems and outbreaks of waterborne diseases [59,60]. Human-induced climate change has increased rainfall associated with tropical cyclones, and the frequency and intensity of extreme precipitation events are projected to rise further [61,62]. As extreme rainfall events become more frequent, human habitats face escalating challenges and threats. Therefore, in addition to general evaluation metrics, this study also explores the improved model’s ability to represent extreme precipitation events [5,12]. Reanalysis systems often fall short in capturing extreme precipitation events due to their coarse spatial resolution and the complexity of parameterizing convective weather systems, which are challenging to accurately represent. Although ERA5, with a resolution comparable to satellite products, outperforms MERRA2 and CFSR during summer in this study, its performance still falls short of being sufficiently reliable for hydrological research during this season. Although the ATW method proposed in this study demonstrated the highest performance among all products during the summer, its performance in this season remains the weakest when compared to other seasons.

A quantile mapping method was proposed to correct RPE’s inherent flaw to represent high precipitation values. Results showed that both ATW and TW products demonstrated improved performance, with ATW significantly enhanced after correction using the high-quality APHRODITE dataset. Specifically, ATW showed a considerable improvement over the initial TW when applied for extreme value correction, achieving better results. However, the QM method was found to be ineffective in fully correcting extreme precipitation values. Despite improvements, the underlying limitations of RPE products, particularly their inability to accurately express high-intensity rainfall, remain a challenge. According to Zhao et al. [63], QM is an unconditional method, meaning it does not preserve the relationship between predicted and observed data pairs. Consequently, QM may sometimes mislead the correction process of the original data. While QM can reduce biases to some extent, it does not ensure full reliability. Despite the limitations of QM for extreme precipitation correction, the improvements in ATW over TW suggest that the multi-source data integration framework proposed in this study is indeed effective.

5. Conclusions

The main conclusions of this study are as follows:

(1) Among the three RPEs, MERRA2 shows the lowest uncertainty, although this is closely tied to the grid’s average rainfall values. In this context, calculating the signal-to-noise ratio (SNR) is useful for assessing the strength of each product relative to background noise. SNR-based weighting indicates that, for mainland China, MERRA2 has the highest weight during summer and autumn, while ERA5 has a higher weight in spring and winter. As seen in Figure 4, MERRA2’s higher SNR in the summer is concentrated in mid-to-low latitudes, whereas ERA5′s elevated winter weight stems from its higher SNR in the Song-Liao River Basin.

(2) The TW product, generated by merging the three RPEs using the TCH method, demonstrates improved alignment with observed data. However, without correcting for rainfall occurrence/non-occurrence, TW fails to capture more nuanced precipitation characteristics. After incorporating corrections from the high-density APHRODITE dataset, the ATW product achieves significantly greater accuracy compared to TW and other individual products. The improvement in ATW is consistent across all watersheds and seasons, substantially enhancing the estimation of the $KGE \prime$ compared to original RPE products.

(3) Despite these advancements, limitations remain due to inherent flaws in RPE products, particularly in their ability to represent extreme precipitation events. While the QM method improves the alignment of extreme rainfall values with observations, certain discrepancies persist.

The RPE integration methods and products proposed in this study can be applied to climate change research and hydrological modeling. Additionally, the main findings will be valuable for research aimed at enhancing reanalysis systems. Future research should focus on refining these integration techniques to better capture the nuances of extreme precipitation events and reduce uncertainty. Exploring additional data sources, such as high-resolution observational datasets and next-generation satellite products, can further enhance accuracy. Addressing inherent limitations of current reanalysis products, while integrating innovative statistical and machine learning approaches, will contribute to more reliable precipitation estimation, with broad implications for climate modeling, water resource management, and disaster risk reduction.