Abstract

Temporally consistent high-quality, high-resolution multisensor precipitation reanalysis (MPR) products are needed for a wide range of quantitative climatological and hydroclimatological applications. Therefore, the authors have reengineered the multisensor precipitation estimator (MPE) algorithms of the NWS into the MPR package. Owing to the retrospective nature of the analysis, MPR allows for the utilization of additional rain gauge data, more rigorous automatic quality control, and post factum correction of radar quantitative precipitation estimation (QPE) and optimization of key parameters in multisensor estimation. To evaluate and demonstrate the value of MPR, the authors designed and carried out a set of cross-validation experiments in the pilot domain of North Carolina and South Carolina. The rain gauge data are from the reprocessed Hydrometeorological Automated Data System (HADS) and the daily Cooperative Observer Program (COOP). The radar QPE data are the operationally produced Weather Surveillance Radar-1988 Doppler digital precipitation array (DPA) products. To screen out bad rain gauge data, quality control steps were taken that use rain gauge and radar data. The resulting MPR products are compared with the stage IV product on a daily scale at the withheld COOP gauge locations. This paper describes the data, the MPR procedure, and the validation experiments, and it summarizes the findings.

1. Introduction

Accurate historical precipitation analysis is needed for various hydrologic, hydrometeorological, and hydroclimatological applications. Increasingly, many of these applications require analysis at high (4 km and hourly or higher) spatiotemporal resolutions. Since the implementation in the early to mid-1990s, the network of Weather Surveillance Radar-1988 Doppler (WSR-88D), commonly known as the Next Generation Weather Radar (NEXRAD), real-time precipitation analysis in the United States has been heavily relying on radar data for high-resolution precipitation information. In the continental United States, the WSR-88D network consists of approximately 140 sites, most of which have been operational for over a decade. The WSR-88Ds provide radar reflectivity estimates for the NEXRAD precipitation processing subsystem (PPS; Fulton et al. 1998), which produces radar-derived precipitation products in real time in support of the National Weather Service’s (NWS) mission and external users. Quantitative precipitation estimates (QPE) from radar, however, are subject to various sources of error (Wilson and Brandes 1979; Vasiloff et al. 2007) and, by themselves, are generally not suitable for quantitative hydroclimatological applications. To produce precipitation estimates that are more accurate than those obtainable from radar or rain gauges alone, multisensor estimation is necessary, consisting usually of bias correction of radar QPE and multivariate analysis, or merging, of bias-corrected radar QPE and rain gauge data. In NWS, such multisensor precipitation estimation (MPE) applications produce precipitation estimates at different spatial scales and in stages (Hudlow 1988; Vasiloff et al. 2007). For example, the so-called stage III products are generated at the River Forecast Centers (RFCs), which are nationally mosaicked at the National Centers for Environmental Prediction (NCEP) to produce the stage IV products (Fulton et al. 1998; Vasiloff et al. 2007). In the early 2000s, the MPE algorithm replaced the stage III algorithm at most RFCs (Breidenbach et al. 1998, 2001a, 2002). At the Arkansas–Red Basin River Forecast Center and in the western region, the Process 1 algorithm (see appendix A of Seo and Breidenbach 2002) and the Mountain Mapper (see section 2 of Schaake et al. 2004) are used, respectively, instead of the MPE.

Since they were first implemented, both the radar and multisensor QPE algorithms in PPS and MPE, respectively, have undergone a number of improvements, resulting in significant changes in the error characteristics of the respective precipitation products. Being a real-time operation, stage III analysis can only make use of the real-time data that are available by the time of the analysis. The objective of this work is to produce temporally consistent high-quality, high-resolution multisensor precipitation reanalysis (MPR) products for a wide range of climatological and hydroclimatological applications. Development of such products should capitalize on the additional data that may not be available in real time and on the retrospective nature of the analysis that allows for identification, correction, and quantification of the errors, which are generally not possible in real time due to insufficient data and computing power. Therefore, the National Climatic Data Center (NCDC), in collaboration with the NWS, has redeveloped the MPE algorithm (Seo and Breidenbach 2002) into the MPR algorithm over the pilot domain (Fig. 1). MPR inputs the historical digital precipitation array (DPA) products from the WSR-88D network and the rain gauge data from the Cooperative Observer Program (COOP) and the Hydrometeorological Automated Data System (HADS) networks and outputs a suite of reanalysis products similar to that of MPE products (see Fig. 2). As a pilot project, the reanalysis is set up for the regional domain over North Carolina and South Carolina, which includes six WSR-88D sites (see Fig. 1, Table 1). The goal of the pilot project is to demonstrate the improvement of experimental MPR products over the operational QPE products. The main sources of improvement include additional rain gauge data, systematic quality control (QC) of rain gauge data, correction of systematic biases in radar QPE, and parameter optimization for radar–rain gauge merging. In this paper, we describe the data and the reanalysis procedure used for the pilot project and summarize the results, including comparative evaluation with the stage IV products. This paper is organized as follows: section 2 describes the data sources used in the study; section 3 describes the MPR procedure; section 4 provides an overview of the cross validation of the MPR procedure; and section 5 provides conclusions and future research recommendations.

Fig. 1.

Study location including radar range rings of 230 km and rain gauge locations.

Fig. 1.

Study location including radar range rings of 230 km and rain gauge locations.

Fig. 2.

Workflow of the MPR.

Fig. 2.

Workflow of the MPR.

Table 1.

WSR-88D sites used in study.

WSR-88D sites used in study.
WSR-88D sites used in study.

2. Data and quality control

The radar data source is the WSR-88D DPA and the rain gauge data sources are HADS and COOP. In this section we provide an overview of the data and the quality control steps taken to screen out bad data and gauges.

a. DPA

The DPA is 1 of approximately 20 products in the WSR-88D level III data suite that are derived from the level II data (Klazura and Imy 1993). The DPA product is an hourly running precipitation total on the Hydrologic Rainfall Analysis Project (HRAP) grid (approximately 4 × 4 km2) out to a range of 230 km from the radar. The estimates are produced by the WSR-88D PPS (Fulton et al. 1998) by combining radar reflectivity estimates from multiple elevation angles using a hybrid scan strategy, converting reflectivity Z to rainfall R using the ZR relation and then mapping them onto the HRAP grid. Currently, five sets of parameters are used operationally for the ZR relation (Z = ARb) for five different precipitation regimes: A = 300, b = 1.4; A = 200, b = 1.6; A = 250, b = 1.2; A = 130, b = 2.5; and A = 75, b = 2.5 for convective, stratiform, tropical, stratiform-east, and stratiform-west precipitation events, respectively. The major sources of errors in the DPA include, but are not limited to, radar calibration (Anagnostou et al. 2001; Ulbrich and Lee 1999; Smith et al. 1996), anomalous propagation (Krajewski and Vignal 2001), brightband enhancement (Smith et al. 1996), radar beam blockage (Young et al. 1999), nonuniform vertical profile of reflectivity (Seo et al. 2000), and uncertain microphysical parameters, such as the ZR coefficients (Smith et al. 1996; Vasiloff et al. 2007). These errors vary in space and time and accumulate nonlinearly in time. As such, it is extremely difficult, if not practically impossible, to isolate individual errors from the DPA products and correct them a posteriori. In this work, we make no attempt to improve the intrinsic quality of the DPA products as, in our view, such an effort is more cost effective by reprocessing the level II data. As such, we rely solely on the operationally produced DPA products for radar QPE.

To screen out obviously bad DPA data from the reanalysis process, we developed an ad hoc procedure. It consisted of first calculating the fractional coverage and conditional mean of radar precipitation for each of the hourly rainfall maps associated with each of the six radars in the pilot domain. Next, we defined a subjective threshold value for each of the statistics and identified the hours for which the statistics exceeded this threshold. We then created an animation of these hours and visually identified the obviously bad DPA precipitation maps that could be attributed to ground returns from anomalous propagation, bright band, or other nonprecipitating events, such as birds and insects. These hours were then added to a list that was used as a database for exclusion during the follow-on processing of the MPR algorithm. While this ad hoc procedure is possible for the regional domain, we recognize that a procedure like this may not be fully implementable at the conterminous United States scale. In addition, procedures for automated QC of radar data are relatively immature, even at the radar operations level. At this time we believe this is another challenge that is identified by the reanalysis work, which we identify as a project that could stand on its own.

b. HADS

The primary source of hourly rain gauge data are the reprocessed HADS data of Kim et al. (2009), which are based on the historical archive of the real-time data collected and distributed by the HADS Program in the NWS Office of Hydrologic Development (OHD). The period of record of the reprocessed data is 2002–present. To the best of the authors’ knowledge, HADS contains the only hourly precipitation dataset available for the entire United States as a single product. Kim et al. (2009) have investigated data loss and quality issues with HADS, including possible biases with respect to neighboring COOP gauges. Through reprocessing, they addressed a number of these issues, which resulted in an increase in the number of hourly HADS precipitation data (particularly those ending at the top of the hour) and improved quality. In many instances, the reprocessed data resulted in the recovery of storm events and other hourly observations that may be critical to MPR. For further details, the reader is referred to Kim et al. (2009), who present many examples of data quality issues and some stumbling blocks associated with the HADS dataset.

In addition to the rain gauge database quality control steps used in reprocessing HADS data (see Kim et al. 2009 for details), we have implemented a simple procedure to identify bad rain gauges using radar data. We screened out bad gauges by comparing seasonal accumulations of rain gauge precipitation and the corresponding radar pixel values in the DPA products. We removed suspect rain gauges by visual inspection and examined a set of summary statistics for individual gauges for the given season. The statistics include the indicator and conditional correlation coefficients to check the strength of association in precipitation detection and estimation of precipitation amount given detection, respectively, between the rain gauge observation and the collocated radar estimate within the effective coverage of the radar (Breidenbach et al. 1999, 2001b) and the conditional coefficient of variation of precipitation amount to check the integrity of its seasonal statistics against climatology. For visual inspection, we also examined the scatterplots of rain gauge values versus the corresponding radar pixel values of hourly precipitation at the seasonal scale. This additional quality control results in screening out approximately another 20% of the reprocessed HADS data.

c. Cooperative Observer Program

The Cooperative Observer Program provides thousands of meteorological and hydrometeorological observations daily in the United States (NCDC 2009). These data include both hourly and daily gauge precipitation observations that can be used in MPR. The spatial density of the daily stations is much larger than that of the hourly stations. As such, it is important that MPR utilize the daily COOP observations to the maximum extent possible. Because the analysis in MPR is currently carried out only at the hourly scale, it is necessary to disaggregate the daily gauge observations into hourly estimates of gauge precipitation. Section 3c describes the disaggregation procedure used in MPR.

We have applied the DPA-based quality control procedure described in section 2b to identify suspect gauges that may have survived the gauge database quality control steps. As with HADS data, the procedure was carried out for each season of the analysis period. This additional quality control screens out approximately 10% of the COOP data. Table 2 summarizes the availability of HADS and COOP rain gauge data in the pilot domain.

Table 2.

Rain gauge network characteristics in the pilot domain. Density units are square kilometers per gauge.

Rain gauge network characteristics in the pilot domain. Density units are square kilometers per gauge.
Rain gauge network characteristics in the pilot domain. Density units are square kilometers per gauge.

3. Reanalysis procedure

The reanalysis procedure includes several steps that are designed to take full advantage of the retrospective nature of the analysis. They include incorporating additional rain gauge data that were not available in the real-time analysis, systematic quality control of the input data, correction of systematic biases in the radar precipitation data, and optimization of the parameters used for multisensor precipitation estimation. Figure 2 shows the MPR workflow that identifies these steps and the input to and output from them, along with the products generated. Next we describe each of the steps in some detail.

a. Mosaicking of radar QPE

After removing the obviously bad data, the DPA is used to create a mosaic of the radar-only QPE data over the analysis domain, referred to as Rmosaic. The procedure used here is the same as that in MPE (Breidenbach et al. 1999, 2001b), which involves identifying the effective coverage of the radar based on the long-term seasonal radar climatology and the known height of the radar beam. The mosaicking process consists of the following steps: For each HRAP bin in the analysis domain, determine if the bin falls within the effective coverage of any radar. If it does not, the radar QPE is not available for that bin. If it falls within the effective coverage of two or more radars, identify the bin with the lowest unobstructed sampling height and assign the radar precipitation estimate from that radar to that bin. For further details, the reader is referred to Breidenbach et al. (1999, 2001b). While less than optimal, this estimation procedure identifies very well the significant interradar mean field biases that arise from various sources. Next we describe how such biases are corrected in MPR.

b. Correction of mean field bias

In Rmosaic, mean field bias manifests as artificial linear boundaries equidistant from any two adjacent radar locations when precipitation is accumulated over time. Figures 3a and 3b show the precipitation accumulations over a 6-month period for two adjacent radars. Note that, over a season, the differences can be seen clearly. Some of the differences are due to different choices of certain adaptable parameters in the precipitation processing subsystem (e.g., the ZR coefficients), the vertical profile of reflectivity effects, and beam blockage. For long-term accumulations, however, the interradar calibration difference (Anagnostou et al. 2001; Ulbrich and Lee 1999; Smith et al. 1996) is usually the main contributing factor. Figure 3c shows the difference, expressed as the ratio between the two radar precipitation estimates for the same grid boxes. Note that there are areas in the overlapping region with a ratio of 2:1 or greater.

Fig. 3.

WS accumulation (mm) of the radar DPA at (a) Columbia, SC, and (b) Wilmington, NC, and (c) the resulting ratio of Wilmington to Columbia in the overlap area of the two radars.

Fig. 3.

WS accumulation (mm) of the radar DPA at (a) Columbia, SC, and (b) Wilmington, NC, and (c) the resulting ratio of Wilmington to Columbia in the overlap area of the two radars.

We correct for interradar calibration differences by comparing, radar site by radar site, DPA estimates to rain gauge observations in the long term. This correction is conceptually similar to the mean field bias correction in MPE at a large time scale (Seo et al. 1999; Smith and Krajewski 1991) as explained next. We define the mean field bias as the ratio of the long-term accumulation between the COOP data within the effective coverage of the radar (Breidenbach et al. 1999, 2001b) and the collocating DPA estimates from that radar

 
formula

where N denotes the number of COOP gauges within the effective coverage of the radar for that particular site and season. Throughout this paper, we define the warm season (WS) as April–September and the cool season (CS) as October–March. We delineate these seasons as warm and cool, because studies have shown that when the convective season begins, the reflectivity climatology changes almost instantaneously (Steiner 1996; Atlas et al. 1990). Figure 4 shows the mean field bias for each radar and season as estimated via Eq. (1). Notice that the bias has reduced (i.e., closer to unity) over time in large part to the improvements made to the precipitation processing subsystem.

Fig. 4.

Long-term bias values of COOP rain gauge locations at each corresponding radar umbrella for the DPA pixel locations for (a) CS and (b) WS.

Fig. 4.

Long-term bias values of COOP rain gauge locations at each corresponding radar umbrella for the DPA pixel locations for (a) CS and (b) WS.

The mean field bias thus obtained is applied to the DPA estimates for the specified season for each radar site. Then, the adjusted estimates are mosaicked in the same way as Rmosaic to create the bias-adjusted hourly precipitation field Bmosaic over the analysis domain (Fig. 2). Figure 5 shows the result of the mean field bias correction at the seasonal scale. Note that the linear boundaries are now removed.

Fig. 5.

Merged seasonal accumulation (mm) of (a) DPA values and (b) interradar bias-corrected values as described in section 3a for the 2005 WS.

Fig. 5.

Merged seasonal accumulation (mm) of (a) DPA values and (b) interradar bias-corrected values as described in section 3a for the 2005 WS.

c. Time disaggregation of daily COOP

To utilize all available rain gauge data in the pilot domain in hourly analysis, we disaggregate daily rain gauge observations to hourly rain gauge estimates. The procedure used in this work is based on a similar procedure used for MPE in NWS operations (NWS 2009). It involves pairing daily rain gauge observations with collocated DPA estimates. Then, the daily rain gauge estimates are time-distributed proportionally according to the hourly DPA estimates over that 24-h period, which produces hourly estimates of rain gauge precipitation at the COOP locations. Our experience is that using radar precipitation estimates to distribute the COOP data works better for the warm season than for the cool season likely due to larger errors in radar QPE in the cool season. For this reason, we extended the time-disaggregation procedure to use the gauge-only analysis and hourly HADS gauge data to time-distribute the daily COOP observations in the cool season.

Undoubtedly, there are a number of significant questions regarding the accuracy of the resulting hourly estimates of rain gauge precipitation. Indeed, one of the motivations of MPR is to assess as a part of the planned product development and improvement the effect of using such hourly estimates of gauge precipitation on hourly or subdaily precipitation analysis. It is also worth noting that, because the multisensor algorithm used in MPR is an exact interpolator (Seo 1996), the daily analysis based on summing the 24 hourly analyses using the hourly estimates of gauge precipitation reproduces the observed daily gauge amounts at the daily gauge locations. At other locations, however, the 24-h sum of hourly analyses would differ from daily analysis based on 24-h gauge accumulations for the same period because of the nonlinearities involved. In the section 4, we assess the effect of using time-disaggregated hourly estimates of gauge precipitation in hourly analysis for 24-h precipitation.

d. Local bias correction

The Bmosaic field, which is obtained from mosaicking the mean field bias-corrected DPA products from multiple sites as described earlier, may have spatially varying biases that may be correctable by rain gauges, depending on the quality of the radar QPE and the density of the available rain gauges. The sources of such “local” biases include partial beam blockage due to structures and terrain (Young et al. 1999), returns from anomalous propagation of the radar beam, and nonuniform vertical profiles of reflectivity, including brightband enhancement (Baeck and Smith 1998; Smith et al. 1996; Young et al. 1999). The local bias-correction procedure used in MPR is conceptually similar to that used in MPE (Seo and Breidenbach 2002) but applied at a seasonal time scale. The first part of the procedure is to optimize the radius of influence via cross validation using the following steps:

  1. Assume a starting radius of influence to be used to collect collocated radar–rain gauge pairs at the seasonal scale for evaluation of local biases at gauge locations.

  2. Calculate the gauge-to-radar biases at all gauge locations using seasonal accumulations of the collocated radar–rain gauge pairs.

  3. Withhold each rain gauge-to-radar bias and, using the rest of the rain gauge-to-radar bias values within the radius of influence, estimate via spatial interpolation the rain gauge-to-radar bias at the withheld location.

  4. Apply the estimated rain gauge-to-radar bias to hourly Bmosaic at the withheld location and evaluate the objective function (see next).

  5. Repeat steps 4 and 5 for all rain gauge locations and evaluate the objective function.

  6. Repeat steps 1–6 for a range of values for radius of influence (Fig. 6). Choose the smallest radius of influence associated with the minimum objective function value.

Fig. 6.

Objective function (mm; RMSE) vs. radius of influence (km) as described in section 3b for one seasonal analysis as an example of the parameter optimization for the local bias-correction procedure.

Fig. 6.

Objective function (mm; RMSE) vs. radius of influence (km) as described in section 3b for one seasonal analysis as an example of the parameter optimization for the local bias-correction procedure.

In step 4, the objective function consists of RMSE of hourly precipitation estimates and RMSE of the quantile–quantile (QQ) plot of hourly precipitation estimates relative to the one-to-one line. The motivation for the second term is to reduce conditional biases (e.g., under- and overestimation of large and small precipitation amounts, respectively) that are typical in the estimation of noisy and skewed variables (Fuller 1987). Such an objective function has been used successfully in Seo et al. (2006).

Once the radius of influence is estimated from these six steps for each season, local bias correction is performed on the hourly Bmosaic fields for that season as follows:

  1. Accumulate the gauge and collocated Bmosaic values seasonally and calculate the rain gauge-to-radar bias values at all gauge locations. This method follows the same concept of Smith and Krajewski (1991) and Seo and Breidenbach (2002), in that using the multiplicative bias we can cover most of the variations in the local bias.

  2. Spatially interpolate the rain gauge-to-radar bias values using the optimum radius of influence. This procedure results in a local bias map for the season.

  3. Apply the local bias map to the hourly Bmosaic fields from section 3a to obtain hourly local bias-corrected fields (Lmosaic in Fig. 2) for the entire analysis domain for the season.

Figure 6 shows an example of the resulting radius of influence as a function of radius of influence. In this example, the warm- and cool-season minima exist at about 160 and 240 km, respectively, though significant sampling uncertainties may be seen. We also note here that the optimum radius of influence thus estimated shows significant interannual variability.

One may expect that the quality of Lmosaic would be better if local bias could be estimated reliably at a time scale smaller than seasonal, so that one may resolve spatially varying biases at an event or synoptic scale. We note here that we tried smaller time scales for estimation of local bias, ranging from 5 days to 1 month, instead of 1 season. The results indicate, however, that sampling uncertainty is too large for smaller scales and that the seasonal scale provides the best overall performance.

e. Radar–rain gauge merging

The bias-correction steps described earlier are intended to reduce mean error but not necessarily error variance. The primary purpose of radar–rain gauge merging is to reduce error variance. The merging algorithm used for MPR is a variant of the operational procedure used in MPE by the NWS, referred to as single optimal estimation (SOE; Seo 1998). The algorithm has many adaptable parameters. For MPR, we have identified the following four parameters to be optimized for each month: the multiplicative bias factor, the indicator and conditional correlation coefficients between radar and rain gauge precipitation, and the coefficient of variation of precipitation. The first step in the merging process is to optimize these parameters on a monthly basis as follows:

  1. Specify the a priori (i.e., the best guess) parameter values.

  2. Evaluate, hour by hour, via cross validation the merged estimate at each gauge location using a variant of single optimal estimation (Seo 1996).

  3. Evaluate the objective function made of a combination of error statistics (see next) of the merged estimate against the gauge observation.

  4. Evaluate as necessary the gradient of the objective function with respect to the four parameters being optimized.

  5. Repeat steps 2–4 until convergence is reached.

The algorithm used for the optimization is Fletcher–Reeves–Polak–Ribiere minimization (Press et al. 1992), a conjugate gradient method. The gradients were calculated from the tangent linear code generated from Tapenade (available online at http://tapenade.inria.fr:8080/tapenade/index.jsp). Once the parameters are optimized for each month of each year, the merging algorithm is run using the optimized parameter values for the entire reanalysis period.

The merging algorithm also estimates the estimation variance associated with the precipitation estimate (Seo 1996). Owing to the cross-validation procedure used in parameter optimization, it is possible in MPR to evaluate the quality of the estimation variance estimated from single optimal estimation. We only note here that the estimates tend to be biased because of the limited modeling of heteroscedasticity, but he bias can be removed post factum by applying a multiplicative correction factor.

4. Evaluation

In this section we describe how the MPR products are comparatively evaluated among themselves and with the stage IV product (Vasiloff et al. 2007), which is arguably the most accurate operational QPE product today.

a. Cross-validation experiment

To evaluate the MPR products in the Carolinas’ domain, we designed and carried out an experiment to cross validate the MPR estimates at a daily scale at daily COOP locations. The experiment is designed to allow a fair comparison with the stage IV product, as explained later. The stage IV product is a mosaic over the contiguous United States of the stage III products produced by the NWS RFCs. The stage III product over the pilot domain is generated by the Southeast River Forecast Center (SERFC) and mosaicked into the stage IV product at NCEP. The stage III product uses HADS hourly gauge data. As such, comparison of the MPR product with the stage IV at HADS locations would necessarily favor the latter. On the other hand, the stage III product over the pilot domain did not use the daily COOP data. Accordingly, cross validation of the MPR product at the COOP locations against the daily COOP observations provides a fair comparison with the stage IV product at those locations. It also meant, however, that evaluating the MPR estimates at subdaily scales, such as hourly and 6 hourly, was not possible.

The performance measures used in the comparative evaluation are the bias, the correlation, and the mean squared error and its decomposition of the daily estimates as compared with the “truth” (i.e., daily rain gauge measurements). We define the bias as the ratio of the estimate to the truth

 
formula

where Gk is the gauge observation at some COOP gauge location for day k, Rk is the radar estimate at the same location for the same 24-h period, and K is the total number of days in the evaluation period. Performance for estimation of large precipitation amounts is particularly important for many applications. As such, we evaluate and present both the unconditional (≥0 mm) and conditional (>6.4, >12.8, and >25.4 mm) statistics. These thresholds ensure that there is a large enough data sample size for computing statistics.

b. Results

Figures 7 and 8 show the cross-validation results for warm and cool seasons, respectively, in the form of scatter and QQ plots [see Kottegoda and Rosso (1997) for a description of quantile–quantile plot]. The figures show the verifying gauge observation versus the estimate [Gmosaic, MPR product (Mmosaic), or stage IV] for daily amounts for warm or cool seasons for the period 2002–07. The QQ plots show how closely the marginal distribution of the estimated precipitation matches that of the observed. The discrepancies between the gauge and radar estimates (i.e., the scatter) are in part due to the natural variability of precipitation in space and time. The figures show that Mmosaic has smaller scatters than stage IV for both warm and cool seasons but slightly larger scatter than Gmosaic for the cool season. The QQ plots show that Mmosaic most closely follows the marginal distribution of observed precipitation for larger precipitation amounts. The small number of highly overestimated stage IV values in the cool season is probably due to bad HADS data that could not be quality controlled in real time.

Fig. 7.

Scatterplots of product estimates for (a) Mmosaic, (b) stage IV product, (c) Gmosaic, and (d)–(f) corresponding QQ plots for all WS daily scale estimates vs rain gauge estimates from 2002 to 2007. All units are millimeters.

Fig. 7.

Scatterplots of product estimates for (a) Mmosaic, (b) stage IV product, (c) Gmosaic, and (d)–(f) corresponding QQ plots for all WS daily scale estimates vs rain gauge estimates from 2002 to 2007. All units are millimeters.

Fig. 8.

As in Fig. 7, but for CS.

Fig. 8.

As in Fig. 7, but for CS.

Figure 9 shows the bar graph of the bias [Eq. (2)] for each MPR product and the stage IV product. In the figure, the Gmosaic is the gauge-only estimate (Fig. 2). Figure 9 shows that the biases are reduced in each step of the MPR process, from Bmosaic to Lmosaic to Mmosaic. The stage IV product has little bias in the warm season in the pilot domain, but it has a significant bias in the cool season. The major factors contributing to cool-season bias are radar range dependency (Smith et al. 1996; Young et al. 2000), beam overshoot of low-level stratiform events (Young et al. 1999), ZR relation biases due to frozen precipitation (Hunter and Holroyd 2002), and brightband enhancement (Smith et al. 1996). This cool-season bias is due to the large low bias (i.e., underestimation) in Rmosaic, and it reflects the difficulty of completely correcting it in real time using limited rain gauge data. Mmosaic, on the other hand, is essentially bias free for both warm and cool seasons.

Fig. 9.

Long-term bias between the rain gauge values and the algorithm products—Gmosaic, Mmosaic, stage IV product, Lmosaic, Bmosaic, and Rmosaic.

Fig. 9.

Long-term bias between the rain gauge values and the algorithm products—Gmosaic, Mmosaic, stage IV product, Lmosaic, Bmosaic, and Rmosaic.

Figures 10a and 10b show the unconditional and conditional RMSE and correlation for each estimate for warm and cool seasons, respectively, from 2002 to 2007. The smaller the RMSE and the larger the correlation coefficient is, the better the estimate is. The following may be seen in Fig. 10. In the warm season, Mmosaic is consistently superior to stage IV, and stage IV is consistently superior to Gmosaic. For daily amounts greater than 25.4 mm, the marginal improvement of Mmosaic over Gmosaic is, as expected, much greater. In the cool season, Gmosaic is consistently superior to both Mmosaic and stage IV, and Mmosaic is consistently superior to stage IV. Confidence intervals (5% and 95%) for the >25.4-mm threshold for Mmosaic are 0.016 for correlation and 0.54 for RMSE in the warm season and 0.17 for correlation and 0.38 for RMSE in the cool season. The differences at the lower thresholds are also statistically significant. The confidence intervals for the >25.4-mm threshold for Gmosaic (correlation: 0.016 WS and 0.014 CS; RMSE: 0.66 WS and 0.46 CS) and stage IV (correlation: 0.16 WS and 0.22 CS; RMSE: 0.45 WS and 1.38 CS) as well as for all thresholds are statistically significant. Note that both Gmosaic and Mmosaic greatly outperform stage IV. That Mmosaic is inferior to Gmosaic is an indication that the quality of radar QPE (i.e., the DPA products) must be improved considerably in the cool season to provide value to radar–rain gauge estimation. That Mmosaic very significantly improves over stage IV reflects the value of utilizing all available rain gauge data in precipitation analysis, and it reflects the benefits of reanalysis for more effective bias correction and merging.

Fig. 10.

Correlation vs RMSE (mm) for unconditional, gauge > 6.4 mm, gauge > 12.7 mm, and gauge > 25.4 mm for (a) WS and (b) CS at the daily scale.

Fig. 10.

Correlation vs RMSE (mm) for unconditional, gauge > 6.4 mm, gauge > 12.7 mm, and gauge > 25.4 mm for (a) WS and (b) CS at the daily scale.

One of the motivations for MPR is to assess the value of utilizing additional data sources. Here, we provide a comparative evaluation of the MPR products among different choices of rain gauge networks used—that is, HADS only, COOP only, and HADS and COOP combined. Recall that the daily COOP data are used in this evaluation via disaggregation into hourly rain gauge estimates as described in section 3c. Figure 11 shows the correlation and RMSE of daily precipitation for Gmosaic, Mmosaic, and stage IV for each of the three different gauge networks. Note that, because the evaluation is made at a daily scale, the large uncertainty associated with disaggregating daily COOP observations into hourly estimates of gauge precipitation is not as well reflected in Fig. 11 as it would have in the evaluation at an hourly scale. Nevertheless, using the stage IV statistics as reference, it is easy to see the large additional value of using both the HADS hourly observations and the hourly estimates of gauge precipitation from disaggregating the daily COOP observations over using only one of the two networks. It is interesting to note that, in the warm season, the margin of improvement by Mmosaic over stage IV is greater when the HADS network is used over the COOP network, even though the quality of Gmosaic is higher with the COOP network than with HADS. This apparent contradiction may be a reflection of the large uncertainty associated with disaggregating daily COOP observations using the DPA data. In the cool season, Gmosaic using the COOP data (i.e., daily observations disaggregated into hourly estimates using hourly HADS gauge data) provides a larger improvement over stage IV than using HADS data. Recall from section 3 that disaggregation of daily COOP data in the cool season was based on hourly HADS data. As such, the previously-mentioned result probably reflects the fact that the COOP network is denser than the HADS network in the analysis domain (Table 2). Note also that in the cool season, Mmosaic provides improvement over stage IV only when the COOP data are used. That Mmosaic, which uses the reprocessed HADS data, is inferior to stage IV (albeit only marginally), which used the real-time HADS data, suggests that, without interactive quality control by the human forecasters as applied at the RFCs in the generation of the stage III product, the errors in the cool-season DPA products in the pilot domain may be too large to overcome with radar–rain gauge merging alone.

Fig. 11.

As in Fig. 10, but for (a) HADS and COOP inputs, WS; (b) COOP-only input, WS; (c) HADS-only input, WS; (d) HADS and COOP inputs, CS; (e) COOP-only input, CS; and (f) HADS-only input, CS.

Fig. 11.

As in Fig. 10, but for (a) HADS and COOP inputs, WS; (b) COOP-only input, WS; (c) HADS-only input, WS; (d) HADS and COOP inputs, CS; (e) COOP-only input, CS; and (f) HADS-only input, CS.

To further diagnose the relative performance among the MPR and stage IV products, we carried out decomposition of mean square error (MSE; Murphy and Winkler 1987)

 
formula
 
formula
 
formula

where N denotes the total number of data values; fn and on denote the nth estimate and observation, respectively; F and O denote the estimated and observed precipitation; mF and mO denote the mean of the estimated and observed precipitation, respectively; σF and σO denote the standard deviation of the estimate and observation, respectively; var denotes the variance; and ρ denotes the correlation between the estimated and the observed precipitation. In Eq. (3c), the first and second terms measure biases in the mean and in the univariate variability of the estimated precipitation, respectively, and the third term measures the strength of covariation between the estimated and the observed precipitation (smaller values indicate stronger covariation).

The top and bottom plots in Fig. 12 show the decomposition of unconditional MSE according to Eq. (3) for the MPR and stage IV products for the warm and cool seasons, respectively. For the warm season, the bias in the mean [i.e., the first term in Eq. (3c)] is too small to be discernable. Only Gmosaic and Rmosaic show noticeable bias in standard deviation (i.e., in univariate variability). As such, in the warm season, the relative performance under the MSE criterion is determined almost exclusively by the strength of covariation. The figure shows that, under this criterion, Mmosaic clearly performs the best. In the cool season, the picture is very different. Only Rmosaic has sizeable bias in the mean (due to the significant underestimation). The mean field and local bias-correction procedures do reduce this bias in the mean but only at the expense of increased bias in univariate variability and reduced strength in covariation. Mmosaic, Gmosaic, and Stage IV have little bias in the mean and univariate variability. It is clearly seen that Mmosaic substantially improves over stage IV, but it is outperformed by Gmosaic.

Fig. 12.

MSE decomposition for the unconditional case of the algorithm products—Gmosaic, Mmosaic, stage IV product, Lmosaic, Bmosaic, and Rmosaic—for (a) WS and (b) CS for the daily scale estimates and all seasons from 2002 to 2007.

Fig. 12.

MSE decomposition for the unconditional case of the algorithm products—Gmosaic, Mmosaic, stage IV product, Lmosaic, Bmosaic, and Rmosaic—for (a) WS and (b) CS for the daily scale estimates and all seasons from 2002 to 2007.

Figure 13 is analogous to Fig. 12 and shows the MSE conditional on the observed daily precipitation greater than 25.4 mm. In the warm season, all estimates have significant biases in the mean but much smaller biases in the univariate variability. It is interesting to note that Rmosaic, Bmosaic, and Lmosaic have much smaller biases in univariate variability than Mmosaic, Gmosaic, and stage IV, which suggests that radar QPE and bias-corrected radar QPE may capture at least univariate variability of precipitation better than multisensor QPE. It is seen that the bias-correction procedures work as intended in the warm season, in that Bmosaic improves over Rmosaic, and Lmosaic improves over Bmosaic in terms of both univariate variability and covariation. Mmosaic shows comparable biases in the mean and in univariate variability to stage IV, but it noticeably improves the strength of covariation over stage IV. The relative performance in the cool season for estimation of large precipitation amounts is qualitatively similar to the unconditional performance shown in Fig. 12 but with more pronounced contrasts. Note that Rmosaic has a very large bias in the mean, which is removed by Bmosaic and Lmosaic, but at the expense of much larger bias in univariate variability and weakened covariation. Note also that Mmosaic significantly improves over stage IV, but it does not outperform Gmosaic, which again points to the lack of information content in the cool-season radar QPE.

Fig. 13.

As in Fig. 12, but for gauge > 25.4 mm.

Fig. 13.

As in Fig. 12, but for gauge > 25.4 mm.

5. Conclusions

Temporally consistent high-quality, high-resolution multisensor precipitation reanalysis products are needed for a wide range of quantitative climatological and hydroclimatological applications. To address this need, we have reengineered the MPE algorithms of the NWS into the MPR package. Owing to the retrospective nature of the analysis, MPR allows the utilization of additional rain gauge data, more rigorous automatic quality control, and post factum correction of radar QPE and optimization of key parameters in multisensor estimation. As such, one may expect MPR to produce QPE that is more accurate than those obtainable in real time. To evaluate and demonstrate the value of MPR, we designed and carried out a set of cross-validation experiments in the pilot domain of North Carolina and South Carolina. The rain gauge data used in these experiments are from the reprocessed HADS (Kim et al. 2009) and the daily COOP. The radar QPE data are the operationally produced WSR-88D DPA products. To screen out bad rain gauge data, quality control steps were taken that use rain gauge and radar data. The resulting MPR products are compared with the stage IV product at a daily scale at the withheld COOP gauge locations. The stage IV product in the pilot domain is based on real-time HADS and DPA data and interactive quality control of the data and analysis by the forecasters at the NWS SERFC.

The major findings from the experiments are as follows:

  1. Using both the hourly HADS gauge data and the hourly estimates of time-distributed daily COOP data provides very significant improvement in radar–rain gauge and gauge-only analyses of daily precipitation over using only one of the two networks.

  2. In the warm season, Mmosaic is consistently superior to stage IV, and stage IV is consistently superior to Gmosaic. For daily amounts greater than 25.4 mm, the marginal improvement of Mmosaic over Gmosaic is much greater.

  3. In the cool season, Gmosaic is consistently superior to both Mmosaic and stage IV, and both Gmosaic and Mmosaic very significantly outperform stage IV. That Mmosaic is inferior to Gmosaic is an indication that the quality of cool-season radar QPE (i.e., the DPA products) must be improved considerably to provide value to radar–rain gauge estimation. That Mmosaic very significantly improves over stage IV reflects the value of utilizing all available rain gauge data and the benefits of reanalysis for more effective bias correction and multisensor estimation.

The MPR pilot project described in this paper provides a first step toward reanalysis over the United States. While the value of MPR has been demonstrated for the pilot domain, reanalysis at the national scale faces a number of large challenges, many of which are shared also by real-time analysis. To foster a community-wide effort for MPR and to develop and capitalize on the synergism with real-time analysis, we have formed the MPR Working Group. We welcome any comments or suggestions toward realizing MPR at the national scale and participation in the Working Group.

Acknowledgments

The authors thank the NWS/OHD HADS Program for making available the historical HADS rain gauge data used in this work, the NWS/SERFC for very helpful discussions and clarification throughout the course of this work, and the MPR Working Group for critical input at early stages of this work. The second author was supported in part by the NOAA Climate Program Office’s Climate Prediction Program for the Americas Core Project. This support is gratefully acknowledged.

REFERENCES

REFERENCES
Anagnostou
,
E. N.
,
C. A.
Morales
, and
T.
Dinku
,
2001
:
The use of TRMM precipitation radar observations in determining ground radar calibration biases.
J. Atmos. Oceanic Technol.
,
18
,
616
628
.
Atlas
,
D.
,
D.
Rosenfeld
, and
D. B.
Wolf
,
1990
:
Climatologically tuned reflectivity–rainfall rate relations and links to area–time integrals.
J. Appl. Meteor.
,
29
,
1120
1135
.
Baeck
,
M. L.
, and
J. A.
Smith
,
1998
:
Estimation of heavy rainfall by the WSR-88D.
Wea. Forecasting
,
13
,
416
436
.
Breidenbach
,
J. P.
,
D-J.
Seo
, and
R. A.
Fulton
,
1998
:
Stage II and III post processing of NEXRAD precipitation estimates in the modernized Weather Service.
Preprints, 14th Conf. on Interactive Information and Processing Systems (IIPS) for Meteorology, Oceanography, and Hydrology, Phoenix, AZ, Amer. Meteor. Soc., 263–266
.
Breidenbach
,
J. P.
,
D-J.
Seo
,
P.
Tilles
, and
K.
Roy
,
1999
:
Accounting for radar beam blockage patterns in radar-derived precipitation mosaics for River Forecast Centers.
Preprints, 15th Int. Conf. on Interactive Information and Processing Systems (IIPS) for Meteorology, Oceanography, and Hydrology, Dallas, TX, Amer. Meteor. Soc., 5.22. [Available online at http://ams.confex.com/ams/older/99annual/abstracts/1699.htm]
.
Breidenbach
,
J. P.
,
M. A.
Fortune
,
D-J.
Seo
, and
P.
Tilles
,
cited
.
2001a
:
Multi-sensor precipitation estimation for use by River Forecast Centers during heavy rainfall events.
.
Breidenbach
,
J. P.
,
D-J.
Seo
,
P.
Tilles
, and
C.
Pham
,
cited
.
2001b
:
Seasonal variation in multi-radar coverage for WSR-88D precipitation estimation in a mountainous region.
.
Breidenbach
,
J. P.
,
D-J.
Seo
,
P.
Tilles
, and
M.
Fortune
,
2002
:
Multisensor precipitation estimation for use by the National Weather Service River Forecast Centers.
Preprints, 16th Conf. on Hydrology, Orlando, FL, Amer. Meteor. Soc., 2.5. [Available online at http://ams.confex.com/ams/annual2002/techprogram/paper_26937.htm]
.
Fuller
,
W. A.
,
1987
:
Measurement Error Models.
John Wiley and Sons, 440 pp
.
Fulton
,
R. A.
,
J. P.
Breidenbach
,
D-J.
Seo
,
D. A.
Miller
, and
T.
O’Bannon
,
1998
:
The WSR-88D rainfall algorithm.
Wea. Forecasting
,
13
,
377
395
.
Hudlow
,
M. D.
,
1988
:
Technological developments in real-time operational hydrologic forecasting in the United States.
J. Hydrol.
,
102
,
69
92
.
Hunter
,
S.
, and
E. W.
Holroyd
III
,
2002
:
Demonstration of improved operational water resources management through use of better snow water equivalent information.
U.S. Bureau of Reclamation Rep. R-02-02, 75 pp
.
Kim
,
D.
,
B. R.
Nelson
, and
D-J.
Seo
,
2009
:
Characteristics of reprocessed Hydrometeorological Automated Data System (HADS) hourly precipitation data.
Wea. Forecasting
,
24
,
1287
1296
.
Kottegoda
,
N. T.
, and
R.
Rosso
,
1997
:
Statistics, Probability, and Reliability for Civil and Environmental Engineers.
McGraw-Hill, 735 pp
.
Klazura
,
G. E.
, and
D. A.
Imy
,
1993
:
A description of the initial set of analysis products available from the NEXRAD WSR-88D system.
Bull. Amer. Meteor. Soc.
,
74
,
1293
1311
.
Krajewski
,
W. F.
, and
B.
Vignal
,
2001
:
Evaluation of anomalous propagation echo detection in WSR-88D data: A large sample case study.
J. Atmos. Oceanic Technol.
,
18
,
807
814
.
Murphy
,
A. H.
, and
R. L.
Winkler
,
1987
:
A general framework for forecast verification.
Mon. Wea. Rev.
,
115
,
1330
1338
.
NCDC
,
2009
:
Data documentation for data set 3200 (DSI-3200): Surface land daily cooperative summary of the day.
.
NWS
,
cited
.
2009
:
The National Weather Service River Forecast System user manual documentation.
NOAA/NWS/Office of Hydrologic Development Rep. [Available online at http://www.nws.noaa.gov/oh/hrl/nwsrfs/users_manual/htm/formats.php]
.
Press
,
W. H.
,
S. A.
Teukolsky
,
W. T.
Vetterling
, and
B. P.
Flannery
,
1992
:
Numerical Recipes in C: The Art of Scientific Computing.
2nd ed. Cambridge University Press, 994 pp
.
Schaake
,
J.
,
A.
Henkel
, and
S.
Cong
,
2004
:
Application of PRISM climatologies for hydrologic modeling and forecasting in the western United States.
Preprints, 18th Conf. on Hydrology, Seattle, WA, Amer. Meteor. Soc., 5.3. [Available online at http://ams.confex.com/ams/84Annual/techprogram/paper_72159.htm]
.
Seo
,
D-J.
,
1996
:
Nonlinear estimation of spatial distribution of rainfall—An indicator cokriging approach.
Stochastic Hydrol. Hydraul.
,
10
,
127
150
.
Seo
,
D-J.
,
1998
:
Real-time estimation of rainfall fields using rain gage data under fractional coverage conditions.
J. Hydrol.
,
208
,
25
36
.
Seo
,
D-J.
, and
J. P.
Breidenbach
,
2002
:
Real-time correction of spatially nonuniform bias in radar rainfall data using rain gauge measurements.
J. Hydrometeor.
,
3
,
93
111
.
Seo
,
D-J.
,
J. P.
Breidenbach
, and
E. R.
Johnson
,
1999
:
Real-time estimation of mean field bias in radar rainfall data.
J. Hydrol.
,
223
,
131
147
.
Seo
,
D-J.
,
J. P.
Breidenbach
,
R. A.
Fulton
,
D. A.
Miller
, and
T.
O’Bannon
,
2000
:
Real-time adjustment of range-dependent bias in WSR-88D rainfall data due to nonuniform vertical profile of reflectivity.
J. Hydrometeor.
,
1
,
222
240
.
Seo
,
D-J.
,
H. D.
Herr
, and
J. C.
Schaake
,
2006
:
A statistical post processor for accounting of hydrologic uncertainty in short range ensemble streamflow prediction.
Hydrol. Earth Syst. Sci. Discuss.
,
3
,
1987
2035
.
Smith
,
J. A.
, and
W. F.
Krajewski
,
1991
:
Estimation of the mean field bias of radar rainfall estimates.
J. Appl. Meteor.
,
30
,
397
412
.
Smith
,
J. A.
,
D. J.
Seo
,
M. L.
Baeck
, and
M. D.
Hudlow
,
1996
:
An intercomparison study of NEXRAD precipitation estimates.
Water Resour. Res.
,
32
,
2035
2045
.
Steiner
,
M.
,
1996
:
Uncertainty of estimates of monthly areal rainfall for temporally sparse remote observations.
Water Resour. Res.
,
32
,
373
388
.
Ulbrich
,
C. W.
, and
L. G.
Lee
,
1999
:
Rainfall measurement error by WSR-88D radars due to variations in Z–R law parameters and the radar constant.
J. Atmos. Oceanic Technol.
,
16
,
1017
1024
.
Vasiloff
,
S. V.
, and
Coauthors
,
2007
:
Improving QPE and very short term QPF: An initiative for a community-wide integrated approach.
Bull. Amer. Meteor. Soc.
,
88
,
1899
1911
.
Wilson
,
J. W.
, and
E. A.
Brandes
,
1979
:
Radar measurement of rainfall—A summary.
Bull. Amer. Meteor. Soc.
,
60
,
1048
1058
.
Young
,
C. B.
,
B. R.
Nelson
,
A. A.
Bradley
,
J. A.
Smith
,
C. D.
Peters-Lidard
,
A.
Kruger
, and
M. L.
Baeck
,
1999
:
An evaluation of NEXRAD precipitation estimates in complex terrain.
J. Geophys. Res.
,
104
, (
D16
).
19691
19703
.
Young
,
C. B.
,
A. A.
Bradley
,
W. F.
Krajewski
,
A.
Kruger
, and
M. L.
Morrissey
,
2000
:
Evaluating NEXRD multisensor precipitation estimates for operational hydrologic forecasting.
J. Hydrometeor.
,
1
,
241
254
.

Footnotes

Corresponding author address: Brian R. Nelson, NOAA/NESDIS/NCDC, 151 Patton Ave., Asheville, NC 28801. Email: brian.nelson@noaa.gov