## Abstract

Although it is broadly acknowledged that the radar-rainfall (RR) estimates based on the U.S. national network of Weather Surveillance Radar-1988 Doppler (WSR-88D) stations contain a high degree of uncertainty, no methods currently exist to inform users about its quantitative characteristics. The most comprehensive characterization of this uncertainty can be achieved by delivering the products in a probabilistic rather than the traditional deterministic form. The authors are developing a methodology for probabilistic quantitative precipitation estimation (PQPE) based on weather radar data. In this study, they present the central element of this methodology: an empirically based error structure model for the RR products.

The authors apply a product-error-driven (PED) approach to obtain a realistic uncertainty model. It is based on the analyses of six years of data from the Oklahoma City, Oklahoma, WSR-88D radar (KTLX) processed with the Precipitation Processing System algorithm of the NEXRAD system. The modeled functional-statistical relationship between RR estimates and corresponding true rainfall consists of two components: a systematic distortion function and a stochastic factor quantifying remaining random errors. The two components are identified using a nonparametric functional estimation apparatus. The true rainfall is approximated with rain gauge data from the Oklahoma Mesonet and the U.S. Department of Agriculture (USDA) Agricultural Research Service Micronet networks. The RR uncertainty model presented here accounts for different time scales, synoptic regimes, and distances from the radar. In addition, this study marks the first time in which results on RR error correlation in space and time are presented.

## 1. Introduction

As the use of weather radar data in hydrologic applications continues to increase, so does the need to characterize quantitatively the uncertainty of radar-rainfall (RR) estimates. The high level, major sources, and possible practical consequences of this uncertainty have been well recognized and discussed by many researchers (e.g., Wilson and Brandes 1979; Zawadzki 1982; Austin 1987; Fabry et al. 1992; Kitchen and Jackson 1993; Smith et al. 1996; Winchell et al. 1998; Atlas et al. 1999; Steiner et al. 1999; Westrick et al. 1999; Borga 2002; Krajewski and Smith 2002; Krajewski and Ciach 2003; Jordan et al. 2003; Villarini et al. 2007). Nevertheless, the operational RR products still lack any estimate of their uncertainty (Fulton et al. 1998). This leads to a situation in which both those who issue the products and those who use them know that there are often significant errors in the products, but they have no quantitative information about the probable magnitude or statistical structure of these errors. Consequently, there are no mechanisms to account for the uncertainty of RR input in hydrologic forecasting and decision making (Krzysztofowicz 2001). A possible solution to this acute problem is to generate and present the radar products in a probabilistic form that explicitly and comprehensively reflects their uncertainty structure. However, while various studies address RR error propagation via simulation (e.g., Sharif et al. 2004) or focus on selected error sources (e.g., Vignal and Krajewski 2001), our review (Krajewski and Ciach 2003) reveals that there are no general quantitative uncertainty description methods that could be used for this purpose. This study proposes an empirically based mathematical modeling of the RR error structure that can be directly applied to the probabilistic quantitative precipitation estimation (PQPE) based on the operational RR products.

The many sources of uncertainty in RR estimates are frequently enumerated in the above-mentioned literature. In practice, however, we have little or no quantitative data on these separate errors and their interdependencies. Therefore, we pursue the PQPE objective through the product-error-driven (PED) approach. We assume that, for a given RR product, only the combined effect of all error sources can be measured and modeled mathematically in an empirically supported manner. Consequently, the uncertainty modeling in this study consists of thorough analyses of the products based on six years of level II data from the Oklahoma City, Oklahoma, radar (KTLX). The raw radar reflectivity data are processed with the Precipitation Processing System (PPS) algorithm of the Next-Generation Weather Radar (NEXRAD) system (Fulton et al. 1998), and the generated RR products are compared with the ground reference (GR) that approximates the corresponding (concurrent and collocated) true rainfall. In this study, the GR is based on good quality rain gauge data from two networks: the Oklahoma Mesonet and the Agricultural Research Service (ARS) Micronet.

In the past several years, we made a series of efforts to address the problem of RR uncertainty quantification using simplified models and analyzing specialized empirical data. We also designed, deployed, and operated research-quality rain gauge networks in support of these studies. Our work is documented in numerous papers (e.g., Ciach and Krajewski 1999a, b, 2006; Anagnostou et al. 1999; Ciach et al. 2000, 2003; Habib and Krajewski 2002; Krajewski et al. 2003; Ciach 2003; Habib et al. 2004; Gebremichael and Krajewski 2004; Ciach and Gebremichael 2004) and, herein, we build on the experience gleaned from prior research. The challenging task of reliable quantification of the errors in radar-based rainfall estimates consists of two fundamental problems. The first problem pertains to the difficulty in using the point measurements from rain gauges as the GR for evaluating the inherently area-averaged RR products. The second problem concerns comprehensive and, at the same time, concise quantitative description of the discrepancies between radar-estimated and corresponding true rainfall. The later question can be approached through empirically based mathematical modeling of the relationship between these two quantities. Though such a model can be constructed in different ways, in this study, we discuss the version that, we believe, best suits the PQPE objective.

Much of our research has been focused on the first problem stemming from the point-area differences in the GR. Because the spatial sampling areas of the radar and rain gauge sensors differ by at least eight orders of magnitude, rain gauges provide inaccurate approximations of the true areal rainfall corresponding to the RR estimates that we want to evaluate. The simplest way to account for this GR inaccuracy is the error variance separation (EVS) method investigated in Ciach and Krajewski (1999a) and Ciach et al. (2003) and applied to several studies (Anagnostou et al. 1999; Habib and Krajewski 2002; Chumchean et al. 2003). However, this method is limited exclusively to the estimation of error variances and is of little use in the uncertainty modeling necessary for the PQPE application. A much more general approach, the conditional distribution transformation (CDT) method, is investigated in Habib et al. (2004). This method is designed to filter out GR errors to retrieve from the available data the bivariate probability distributions of radar-estimated and corresponding true areal rainfall. Such a distribution is the starting point of the functional-statistical modeling used in the PQPE methodology that we are developing. However, the practical implementation of the CDT method requires specific information on the small-scale conditional spatial correlations in rain gauge–measured rainfall, conditioned on the RR values. As shown in Ciach and Krajewski (2006), such statistics are highly sensitive to the properties of the conditioning scheme. Adequate data required for the CDT implementation can only be acquired with a number of very dense local rain gauge clusters deployed at several distances from the radar. At present, such facilities are unavailable. We also have not yet developed a less demanding method to account for the GR errors. Therefore, in this study, we assume that single rain gauges can provide sufficiently accurate approximations of the areal rainfall averaged over the grids of the operational RR products in NEXRAD (about 4 km by 4 km). This assumption can be justified because, for the time scales considered here (hourly and longer), the spatial rainfall variability in Oklahoma is relatively small (Krajewski et al. 2003; Ciach and Krajewski 2006). It is our judgment that the area-point differences, although not exactly insignificant, do not affect our results in this study in a critical way.

As mentioned above, the central element of our PQPE methodology is a suitable RR error structure model based on a functional-statistical representation of the uncertain relationship between the radar-estimated and corresponding true rainfall. The first of the model’s two components is a deterministic distortion function describing the systematic conditional biases. The second component is a multiplicative stochastic factor quantifying the remaining random errors. The concept of such a functional-statistical representation first appeared in Ciach and Krajewski (1999b) as an idealized analytical model that was later applied to illustrate the possible effects of conditional bias in RR (Ciach et al. 2000). A preliminary empirical uncertainty analysis based on this concept is described in a small sample study performed by Ciach and Gebremichael (2004). In the present paper, we document for the first time the set of comprehensive analyses that are based on a large data sample and are aimed at developing a general and flexible model of RR uncertainty suitable for the PQPE application. The deterministic distortion function and the stochastic factor in the model that we identify here depend explicitly on the values of RR estimates. We estimate them using nonparametric functional estimation tools that do not require consideration of any specific analytical parameterizations prior to the model identification. To assure a degree of generality in our modeling results, we estimate the model components for three seasons, five range zones from the radar, and the averaging time scales from 1 to 24 h. We also compute the spatial and temporal correlation functions of the stochastic component in the RR error model. These characteristics offer insight into the spatiotemporal dependences in the uncertainties and are useful for the simulations of the error process within the PQPE framework.

The paper is divided into seven sections including the introduction. In section 2, we present the main aspects of the modeling methodology that we apply to quantify the RR uncertainty structure. In section 3, we provide a description of the data sample, its quality control, and its stratification into different subsamples. The results of the estimation and parameterization of the RR uncertainty model components are described in section 4. In section 5, we address the question of spatial and temporal dependences in RR uncertainties. Section 6 discusses future extensions of the work and section 7 summarizes the article.

## 2. Modeling methodology

We define the true rainfall as the volume of rainwater that falls on a specified area in a specified interval of time. Since this definition is derived from fundamental physical notions of space and time, the true rainfall is a well-determined and hydrologically meaningful physical quantity. Such a rigorously defined physical quantity is an objective property of the rain process, despite the technical problems with its accurate measurement. Standard operational algorithms of RR estimation, such as the PPS (Fulton et al. 1998), convert the classical radar measurable, that is, the spatial distribution of radar reflectivity factor, into the estimates of hourly rainfall accumulations averaged over specified areal grids. These estimates are particular approximations of the true rainfall over the grid areas, and it is a well-known fact that their accuracy is fairly low (Zawadzki 1982; Austin 1987; Joss and Waldvogel 1990; Fabry et al. 1992; Smith et al. 1996; Westrick et al. 1999; Young et al. 1999; Habib and Krajewski 2002; Krajewski and Smith 2002; Habib et al. 2004). We define all discrepancies between RR values and the corresponding (concurrent and collocated) values of true rainfall as the errors in the RR product.

As a preliminary step in our analysis, we estimate and remove the overall bias in the RR estimates. It is estimated in the standard way using rain gauge data as approximations of the true rainfall:

where *B*_{0} is the bias factor, *R _{g,i}* is the rain accumulated at a rain gauge site over the considered time scale,

*R**

_{r,i}is the corresponding uncorrected RR value, the index “

*i*” denotes a gauge–radar data pair, and the summation is taken over all the data pairs available for the entire sample. The overall bias removal is performed with the following rescaling formula:

where *R _{r,i}* denote the RR values corrected for the overall sample bias. Assuming unbiased measurement error of the rain gauge data, the spatial sampling error of rain gauge rainfall as an approximation of the area-averaged

*R*

_{true}is random (Ciach 2003). Therefore, a large sample of such data allows estimation of

*B*

_{0}with high accuracy. Estimating it and removing the overall bias in an RR product with the rescaling Eq. (2) can be seen as a standard preliminary step in the product uncertainty analysis. Consequently, in the following equations and in the graphs in this paper,

*R*denotes the estimated RR values after removing the overall bias from the data sample, whereas

_{r,i}*R*denotes the corresponding random variable representing statistically the bias-corrected empirical values.

_{r}In quantifying the uncertainties in a given RR product, we assume that they can be completely characterized by the family of bivariate frequency distributions of radar-estimated rainfall, *R _{r}*, and the corresponding true rainfall,

*R*

_{true}, indexed by the major factors that can affect these distributions:

where, respectively, *A* is the spatial averaging scale, *T* is the time scale (accumulation interval), *d* is the distance from the radar, and *S* denotes the synoptic regime of the precipitation system. These distributions can be approximated from large samples of the RR products and the corresponding rain gauge data, provided that the rain gauge measurements are sufficiently accurate representations of *R*_{true}. In practice, this means that the rain gauges were well maintained, their records were thoroughly quality controlled, and the RR errors are considerably larger than the area-point differences for the evaluated RR product.

Once the bivariate distributions defined by (3) are known and described in a tractable mathematical form, we can use them to answer a variety of practically important probabilistic questions. Of particular interest in hydrology are such questions as: “What is the probability that, for given value of *R _{r}* in a specific RR product, the corresponding

*R*

_{true}is greater than some prespecified threshold?” A more general question that, in fact, defines the PQPE problem is, “What are the probability distributions of the likely true rainfall, given the estimated values in the RR product and the other conditions in (3)?” To answer such questions, we construct a model of the

*R*–

_{r}*R*

_{true}relationship below. The model is derived from (3), but has a more specific mathematical structure that allows us to estimate its components based on the available data and parameterize them using parsimonious analytical approximations.

To simplify the notation, let us focus on one resolution, (*A, T*), distance from the radar (*d*), and synoptic regime (*S*), and omit these indices in all the expressions that follow. To describe the *R _{r}*–

*R*

_{true}relationship for any such specific situation, we choose the following functional-statistical representation consisting of two components:

where *R _{r}* and

*R*

_{true}are concurrent and collocated,

*h*(·) is a deterministic distortion function describing the conditional biases depending on

*R*, and

_{r}*e*(·) is a random variable representing the combined effect of all random error sources, the distribution of which can depend on

*R*. In this representation,

_{r}*R*

_{true}is a random variable describing the possible values of the true rainfall that can occur, given an estimated value of the

*R*product. Therefore, the model defined by (4) clearly fits the PQPE objective.

_{r}To fully identify this model based on real data (the sample of the available approximations of the *R _{r}*–

*R*

_{true}pairs), we have to estimate the deterministic distortion function and the statistical distribution of the random component conditioned on

*R*values. Once these components are known, the model defined by (4) can be used to determine the probability distribution of the unknown

_{r}*R*

_{true}for any known

*R*value.

_{r}Because all systematic biases can be described by the overall bias (removed at the beginning of the analysis) and the deterministic distortion function, we can assume without any loss of generality that the expectation of the random uncertainty factor is equal to the unity for each given value of *R _{r}*, that is,

**E**{

*e*(

*R*)

_{r}**|**

*R*=

_{r}*r*} = 1 for any

_{r}*r*value. For this to be true, it is sufficient that we determine

_{r}*h*(·) formally as the following conditional expectation function:

where we use the common statistical convention in which *R _{r}* is a random variable and

*r*is its specified value. Equation (5) is a general form of statistical regression and indicates a straightforward way to estimate the deterministic distortion function,

_{r}*h*(·), based on the available

*R*–

_{r,i}*R*data pairs. It can be estimated using a nonparametric kernel regression framework (e.g., Hardle 1990; Simonoff 1996). The estimator that we use here is the following moving-window weighted averaging:

_{g,i}where *w _{i}* are the weighting factors and

*k*is the “bandwidth” parameter that governs the size of the averaging window centered geometrically on

*r*.

_{r}The specific kernel regression procedure that we apply in this study is discussed in more detail in Ciach (2003). Here, we describe only its basic features. The averaging weights, *w _{i}*, depend on the positions of the

*R*points within the moving window according to a parabolic function called the Epanechnikov kernel (Hardle 1990). The value of the bandwidth,

_{r,i}*k*, controls the trade-off between the smoothing of the data scatter and the preservation the regression curve details. In this study, we applied the same heuristic approach to its selection as in Ciach (2003). We compared visually the regression outcomes for the bandwidth values ranging from 1.2 up to 2.0 and chose the value of

*k*= 1.5 as a satisfactory compromise between oversmoothing and undersmoothing of the estimated regression curves.

Once the overall bias and the deterministic distortion function are known, we can determine the multiplicative random component as

Although its conditional mean is equal to unity for each value of *R _{r}* =

*r*, its distribution can depend on

_{r}*r*. The first step in identifying this dependence is to estimate the conditional variance of the random factor,

_{r}*σ*

^{2}

_{e}(

*r*) = Var{

_{r}*e*(

*R*)

_{r}**|**

*R*=

_{r}*r*}, as a function of

_{r}*r*. This can be accomplished similarly to estimating the

_{r}*h*(

*r*) function. First, notice that the conditional variance can be defined formally as the following conditional mean:

_{r}which allows us to apply to its estimation the same nonparametric technique that was used in (6) to estimate the deterministic distortion component. Consequently, the estimator is the following moving-window weighted averaging:

where the weights, *w _{i}*, the window-size parameter,

*k*, and the parabolic shape of the averaging mask are consistently the same as those applied to (6).

The conditional mean and the conditional variance function of *e*(*R _{r}*) do not provide information on the shape of the conditional distribution of the random component. To know this in more detail, we also estimate the conditional quantiles of

*e*(

*R*) corresponding to predefined levels of the probability of nonexceedence,

_{r}*p*. These quantiles,

*q*, can be defined formally through the following condition:

_{p}where Pr{·|·} is the operator of conditional probability and *q _{p}* is a function of

*r*for each selected value of

_{r}*p*. To be consistent with the nonparametric scheme used to estimate the

*h*(

*r*) and

_{r}*σ*(

_{e}*r*) functions, we compute the conditional empirical quantiles

_{r}*q*(

_{p}*r*) for each selected value of the

_{r}*p*level and each averaging window (determined by the

*r*variable and fixed

_{r}*k*parameter) using a “weighted-point-counting” procedure. It can be defined mathematically as the following formula:

where *e*(*R _{r,i}*) values are determined by (7), and the value of

*x*that fulfills the condition in the braces (curly brackets) above is determined with the precision of 0.001 by simple search. As before, the weights, the window-size parameter, and the parabolic averaging mask are consistently the same as in the estimation procedures applied to (6) and (9).

For the increasing values of the RR estimates, the size of the averaging window in the nonparametric estimators given by (6), (9), and (11) also increases linearly to compensate partially for the dropping frequency of higher rainfall values. Nevertheless, when we get closer to the highest *R _{r}* values, the sample size quickly decreases, which rapidly reduces the statistical significance of the nonparametric estimates of our model components. For this reason, we limit the range of the weighted averaging in the three estimators above to the values of the conditioning

*R*for which the number of points in the averaging window is no fewer than 100.

_{r}We tested the nonparametric estimation schemes described above using different simulation models with their general structure defined by (4) and they proved quite accurate, provided the data sample is large enough. Based on these experiments, we assessed that the size of the 6-yr-long data sample that we use in this study is large enough to obtain reasonably accurate estimates of the empirical model components in three seasons and five distance ranges. However, a comprehensive quantitative assessment of the robustness and stability of these estimates requires extensive analyses and is beyond the scope of this study.

After estimating the deterministic distortion function and the characteristics of the random factor (conditional variance and a number of quantiles) in the nonparametric form, we can approximate the RR uncertainty structure in a more concise (albeit less accurate) manner using simple parametric models. Because the nonparametric estimates reveal the shapes of the empirical functions in much detail, selecting suitable analytical approximations for them is not difficult. Since the extreme rainfall events are the most important in hydrological practice, special attention is required to reproduce the distribution tails of *e*(*R _{r}*) with reasonable accuracy. Some models that have distinctly different tails are the gamma, lognormal, and beta distributions, and each can lead to different decisions in the PQPE applications based on our results. These aspects of our modeling effort and the obtained results are described in detail in section 4.

Apart from the probability distribution, another important characteristic of the random error is its spatiotemporal dependency structure. In the present study, we were able to address this aspect in a limited scope by estimating the spatial and temporal correlation functions of *e*(*R _{r}*). These results are presented in section 5.

## 3. Data sample and its partitioning

The key element of our empirically supported model development is extensive exploratory data analysis based on large samples of RR products and corresponding sets of reliable and independent rain gauge data. This study is based on six years of such data covering the period from January 1998 to December 2003 that was taken from the observation umbrella of the Oklahoma City Weather Surveillance Radar-1988 Doppler (WSR-88D) radar station (KTLX) in the central part of the state of Oklahoma. Data from the KTLX radar are complemented with surface rainfall measurements from two well-maintained rain gauge networks: the Oklahoma Mesonet and the ARS Micronet. Figure 1 shows the location of the rain gauges and of the KTLX radar. In this part of Oklahoma, midlatitude convective systems dominate the rainfall regime (e.g., Houze et al. 1990; Houze 1993), and the average annual precipitation amount is about 800 mm (32 in.).

The official NEXRAD operational RR estimates that are available over the Internet are produced at the River Forecast Centers (RFC) based on the PPS products generated at the WSR-88D stations. The parameters of the PPS processing are often changed by the station operators according to their judgment of what setting is best in the current situation, but the information about the PPS setting used is not archived. In the RFCs, the raw PPS products undergo a number of corrections subject to the judgment of different forecasters. Based on their experience, they use their visual assessment of the products as well as various additional sources of information such as the available rain gauge data. In addition to the differences generated by forecasters, these archived official RR products have been generated over the years with different versions of the PPS algorithms. Therefore, to assure that our large research sample contains a consistent RR product generated under well-specified conditions and that the corresponding rain gauge data can be used as an independent GR, we began our analyses on the lowest level of the not-processed WSR-88D observations in polar coordinates.

We acquired six years of the level II data collected by the KTLX station from the archives of the National Climatic Data Center (NCDC). We processed all these data through the Build 4 of the Open Radar Product Generator (ORPG) version of the PPS software system, which can mimic the workings of the WSR-88D algorithms used operationally in 2005. We then generated the Digital Precipitation Arrays (DPA), which are 1-h accumulation values averaged over the areas of the Hydrologic Rainfall Analysis Project (HRAP) grid (Reed and Maidment 1999; Fulton et al. 1998). To convert reflectivity into rainfall intensity, we used the PPS default settings (Fulton et al. 1998): the standard NEXRAD *Z*–*R* relationship (*Z* = 300*R*^{1.4}), the value of 53 dB*Z* as the hail threshold (“hail-cup”), the value of 0 mm h^{−1} as the minimum signifying precipitation rate, and no range correction. This step of data processing was very time consuming, and it took over three months to make one pass of all the radar files through the PPS software system. However, this was the only way to obtain RR products that, unlike the archived operational products, are generated under known conditions and are independent of the corresponding rain gauge data.

For the GR approximating the true rainfall corresponding to RR values, we used the measurements from available rain gauges. Two well-maintained networks supplied the GR data: the Oklahoma Mesonet (e.g., Elliott et al. 1994; Brock et al. 1995) and the U.S. Department of Agriculture (USDA) ARS Micronet (e.g., Ciach et al. 2003; Habib et al. 2004). The Oklahoma Mesonet was deployed in 1994 and has more than 110 weather stations distributed almost uniformly throughout the state of Oklahoma. At each station, a single tipping-bucket rain gauge manufactured by MetOne Inc. (e.g., Ciach 2003; Ciach and Krajewski 2006) is installed. The intergauge distances in this network are around 40 km, and the network covers a full range of distances from the KTLX site (from about 25 km to the edge of the radar umbrella at 230 km). The experimental USDA ARS Micronet is a relatively dense network of 42 small weather stations covering the Little Washita watershed that has the area of about 900 km^{2}. In the ARS Micronet, the intergauge distances are approximately 5 km, and each station is equipped with a single rain gauge of the same type as in the Oklahoma Mesonet. The stations are located between 70 and 110 km southwest from the KTLX radar site.

The temporal resolution for rainfall measurements in both of the networks is 5 min, and the 5-min rainfall values are determined based on the number of tips generated by a tipping-bucket rain gauge within this interval. As shown in Ciach (2003), the simple method of converting the tips into rainfall limits the accuracy of rainfall measurements considerably, especially at low rainfall intensities and short time scales. In both networks, the tip counting is performed by the dataloggers at the stations, and more detailed information on the tip times that could be used to improve the accuracy (Ciach 2003) is lost irretrievably. To synchronize the surface GR measurements with the analyzed RR product, we aggregated the 5-min rainfall records into the hourly intervals that match the averaging intervals of the RR estimates. For the hourly and longer time scales that we consider in this study, the tipping-bucket time-sampling errors are not significant.

Apart from the availability of an extensive dataset, the critical issue while using rain gauge data as the GR is their effective quality control (QC; e.g., Steiner et al. 1999; Ciach and Krajewski 2006). This has to be performed in a thorough and inquisitive manner, despite the time and effort that it consumes. The standard QC for all the Oklahoma Mesonet data (Shafer et al. 2000) is performed by the Oklahoma Climatological Survey (OCS) as part of the network operation. However, as far as the ARS Micronet data are concerned, the dataset that we were provided with (via Office of Hydrologic Development of the National Weather Service) did not contain any QC flags. Therefore, we developed and implemented an extensive QC procedure following the methods and flagging system similar to that used by the OCS, but adding to it several new tests. After developing our own QC protocol based on the ARS Micronet data (see Krajewski et al. 2005 for details), we applied it to the Oklahoma Mesonet data as well. To illustrate the seriousness of our effort and its outcome, note that it led the OCS to reprocess the QC for their entire database (C. Fiebrich, OCS, 2005, personal communication). This, of course, does not mean that after our QC the rain gauge datasets are perfectly error free. However, we managed to identify and eliminate many problems with these data that had not been previously detected.

The dependence of the RR uncertainties on the synoptic situation or type of the precipitation system noted in (3) is reflected in their strong seasonal dependence. To investigate this effect, we have partitioned the dataset into three seasons: cold (January, February, March, November, and December), warm (April, May, and October), and hot (June, July, August, and September). As guidance for this specific seasonal division, we used the index of the daily average rain-weighted temperature (*T*_{rw}) defined in the following way:

where *R _{g,i}* and

*T*are the 5-min values of rainfall and temperature measured at the stations, and the summations are over all stations and 5-min intervals in the 24-h periods. Using such an index instead of just the surface temperatures is justified by the fact that the directly averaged temperatures are strongly biased by the majority (about 95%) of the 5-min intervals without any precipitation. We use the seasonal partitioning based on (12) to focus on the rainy situations. In Fig. 2, we show the frequency distributions of the

_{i}*T*

_{rw}index for each of the three seasons. The distribution for the cold season is positively skewed; for the warm one it is fairly symmetric; and for the hot season it is negatively skewed.

It is widely acknowledged that the distance from a radar site has a strong effect on the errors in RR estimates. It can result in severe overestimation due to the presence of a bright band (e.g., Smith 1986; Joss and Waldvogel 1990; Fabry and Zawadzki 1995; Smith et al. 1996) or underestimation at far ranges due to partial beam filling and overshooting of the cloud system (e.g., Kitchen and Jackson 1993; Smith et al. 1996; Young et al. 1999). To capture the range effect in our uncertainty model, we divided the radar umbrella into five range zones (Fig. 1). In Table 1, we summarize their characteristics. In selecting the zones, we tried to even out the number of rain gauges in each zone. However, as is evident in Table 1, zone II includes many more gauges than the other zones. This is due to the location of the ARS Micronet network in this zone. Between each zone, we introduced a 5-km-overlap area in order to avoid erratic transition behavior of various statistics.

The DPA maps of hourly RR accumulations are created by the PPS every time a new volume scan becomes available, that is, approximately every 5–6 min. However, because they represent running hourly accumulations, there is strong statistical dependence between the time-adjacent maps separated by short intervals. Therefore, including all the generated DPA products in the analyses is unnecessary and computationally inefficient. On the other hand, considering only the scans separated by a full hour interval results in a drastic reduction of the sample size. In addition, it causes a loss of information because the maps generated between the full hours can be significantly different from the full-hour maps. Including at least part of them effectively increases the stability of our estimates. As a compromise, we included in our estimation procedures only those DPA maps (and the corresponding rain gauge records) that are separated by at least 20% of the considered temporal resolution. This strategy provides more independent information without compromising our sample size. It applies to estimating the overall bias factor as well as to all the nonparametric functions defined in section 2.

## 4. Estimates of the model components

### a. Overall bias factor B_{0}

For each distance zone and season, we computed the overall bias factor, *B*_{0}, using the estimator defined by (1). Table 2 presents the results for hourly scale only because, for an ideal gap-free data sample, the overall bias does not depend on the averaging interval. In practice, due to the gaps in the data, we can observe small variations of the *B*_{0} estimates for longer time scales, but these differences are negligible. The results in Table 2 show that most of the estimated *B*_{0} factors are smaller than 1, which suggests that the radar systematically overestimates the rainfall with respect to the rain gauges. This overestimation tends to maximize in zones III and IV for the hot season, in zone III for the warm season, and in zones II and III for the cold season, which suggests that it can be attributed to the interception of the bright band. On the other hand, the increase of the *B*_{0} factor over the value of 1 indicates the underestimation in RR products. This situation is most pronounced at far distances during the cold season when events are dominated by shallow, low-level precipitation. Thus, it is most likely due to partial beam filling and overshooting of the cloud systems.

### b. Deterministic distortion function h(·)

The overall bias is a single general measure and does not account for the common case when the overestimation and underestimation can occur in the same situation depending on the RR magnitude. This behavior is called a conditional bias, and it is an important characteristic of the RR estimates (Ciach et al. 2000) as well as all other uncertain hydrometeorological products (Katz and Murphy 1997). For this reason, in our RR uncertainty model defined by (4), the overall bias factor is complemented by the deterministic distortion function, *h*(*r _{r}*), quantifying the systematic effect of conditional bias in the form of the dependence of the true rainfall expected value conditioned on the RR values. We estimated this function using the nonparametric estimator given by (6). The estimates were obtained for four time scales (1, 3, 6, and 24 h), for each of the five distance zones and the three seasons.

In Fig. 3, we show the plots of the nonparametric estimates of the deterministic distortion function, *h*(*r _{r}*), for the 1-h accumulation time scale. For the sake of brevity, we present these results exclusively for the hourly scale; however, the same general shape of the functions also holds for the other three time scales (3, 6, and 24 h; see Krajewski et al. 2005). We compare them below in a concise manner after fitting a simple analytical model to the nonparametric estimates. Figure 3 shows that the

*h*(

*r*) curves tend to bend toward the

_{r}*r*axis for higher RR values. Since the overall bias is removed in each case, this bend has to be compensated for in the region of the very small

_{r}*R*values that are abundant in the data sample but not visible in the plot. This behavior of the deterministic distortion function is especially pronounced for the cold season, most likely due to the presence of more small RR values in the cold season than in the other two seasons. For the distances up to 180 km from the radar (below zone V) and for the warm and hot seasons, the conditional biases in Fig. 3 do not show any significant range dependence, and the curves for zones II–IV have no systematic arrangement. This suggests that the description of the variations in the systematic distortions function can probably be simplified, at least outside the cold season and for distances relevant in quantitative hydrological applications of RR products. One outlying exception is the curve for zone I in the hot season. However, the zone I curve also differs in its behavior from the other zones in the cold season. Most likely, this is an indicator that the sample size for zone I is insufficient to obtain stable nonparametric estimates of the

_{r}*h*(

*r*) function (it has the smallest number of rain gauges).

_{r}All the nonparametric empirical estimates of the deterministic distortion function, *h*(*r _{r}*), that we have obtained in this study can be approximated parametrically with satisfactory accuracy using the following power-law function:

The parameters *a _{h}* and

*b*of this analytical model are adjusted to each nonparametric estimate by minimizing the mean-square error (MSE) between the empirical curve and the analytical model defined above. The results of this parameter estimation are summarized in Table 3 for the four temporal resolutions, the three seasons, and the five distance zones. As expected, the estimates of the exponent,

_{h}*b*, are the smallest for the cold season, for which the bends of the curves in Fig. 3 are the strongest. Also, there is no regular behavior of the parameters in respect to the time scales and distance zones, and the warm and hot seasons yield fairly similar outcomes.

_{h}### c. Standard deviation of random component e(·)

While the overall bias and the deterministic distortion function describe the systematic effects in the *R _{r}*–

*R*

_{true}relationship modeled by (4), the random component,

*e*(

*R*), is a stochastic process accounting for the remaining random uncertainties. Its standard deviation,

_{r}*σ*(

_{e}*r*), is a function of the RR values. We estimated this function for different seasons, distance zones, and time scales using the nonparametric regression estimator given by (9). The estimation results for the 1-h accumulation interval are presented in Fig. 4. They show that, for the warm and hot seasons,

_{r}*σ*(

_{e}*r*) level behaves as can be expected: it tends to be larger as distance increases from the radar site. An exception from this regular behavior is the cold season, where a clear distance pattern is not easily distinguishable.

_{r}In general, all the nonparametric estimates of the *σ _{e}*(

*r*) function that we obtained exhibit the hyperbolic behavior:

_{r}where the asymptotic values of *σ*_{0}* _{e}* depends on the zone, season, and temporal scale. We have parameterized the standard deviation of the random component using the following analytical model:

The parameters *σ _{0e}*,

*a*, and

_{e}*b*are fitted to each nonparametric estimate in the same way as for the

_{e}*h*(

*r*) function in (13). We summarize the obtained parameter estimation results in Table 4 for the different seasons, distances, and temporal resolutions. Comparing the results from the analytical fitting of

_{r}*σ*(

_{e}*r*) in Table 4, the parameters in the approximation (15) depend on the temporal scale and the season; however, no clear regularity in these dependences is visible.

_{r}### d. Conditional quantiles of random component e(·)

We already know that the expected value of the random component is equal to 1 and is independent of RR values, whereas its standard deviation is a function of *r _{r}* that can be approximated parametrically using (15). To characterize the probability distribution of the random component in more detail, we estimated selected quantiles as functions of

*r*using the nonparametric estimator defined by (11). In Figs. 5 –8, we present results of this

_{r}*q*(

_{p}*r*) estimation for the probabilities of nonexceedence equal to

_{r}*p*= 0.90, 0.75, 0.50, 0.25, and 0.10 for different seasons and distance zones. These empirical quantiles are compared in the figures with their analytical approximation described below.

From a visual assessment of the quantile plots, we notice that the medians (0.50 quantiles) tend to be around the value of 1, at least for the larger RR values. This suggests a possible symmetric distribution of the random component. In addition, the pairs of 0.25–0.75 and 0.10–0.90 quantiles are approximately symmetric with respect to the medians. For these reasons, we checked two analytical approximations of the empirical quantiles using two symmetric distribution models: the Gaussian distribution with a mean equal to 1 and standard deviation equal to *σ _{e}*(

*r*) and the two-parameter logistic distribution with the same moments [in the applications, the distributions have to be truncated at the value of 0 so that the negative values of

_{r}*e*(

*r*) are not allowed]. In Figs. 5 –8, only comparisons with the Gaussian case is presented because it better fits the nonparametric curves. The figures show a fairly good agreement between the fitted Gaussian quantiles and the empirical nonparametric estimates (with the cold season being the only possible exception), especially for larger RR values. In Fig. 5, in zones II and III, the 0.10 quantile is equal to 0 for values of RR larger than 10 mm. We investigated this feature and concluded that it could be attributed to undetected (i.e., not removed by the PPS) anomalous propagation returns.

_{r}The fact that Gaussian distribution describes the random error factor quite well should not come as a surprise given the overall structure of our uncertainty model defined in (4). Since we consider only conditional errors, the effects of skewness in rainfall itself are removed. Furthermore, the effects of biases are removed by removing the overall and conditional biases prior to the conditional analysis of the random component. This is in contrast with previous studies that assumed a highly skewed lognormal distribution of the RR errors (e.g., Krajewski and Georgakakos 1985; Smith and Krajewski 1993; Ciach and Krajewski 1999b) and used unrealistic models of the uncertainty structure with little empirical support (see Petersen-Øverleir 2005 for additional discussion).

## 5. Spatiotemporal dependences in the random component

In the previous section, we studied the distributional properties of the stochastic component in our RR uncertainty model, *e*(*R _{r}*), conditioned on given RR values. However,

*e*(

*R*) is a random process that can also vary in space and time. Spatiotemporal dependences in the stochastic component are important characteristics of this random process. They are also necessary for simulating realistic RR uncertainty processes in the error propagation models and for generating the likely true rainfall fields given the observed RR within the PQPE framework. The temporal correlations in RR errors are practically unknown. The few studies that tried to account for them in error propagation analyses (e.g., Hossain et al. 2004) are based on strong model assumptions rather than on empirical evidence. The same lack of knowledge applies also to the spatial correlations. The few existing studies that discuss the possible effects of spatial correlations in RR uncertainties (e.g., Nijssen and Lettenmaier 2004; Villarini et al. 2007) use unverified hypothetical models. In this section, we present the first empirical results of these completely unexplored problems.

_{r}To estimate the temporal and spatial correlation functions of the stochastic component in (4), we use standard Pearson’s correlation coefficient:

where *X* and *Y* are either simultaneous values of *e*(*R _{r}*) in two different points in space or two values of

*e*(

*R*) at the same point separated by different time lags. Because of the close-to-Gaussian nature of the random component, Pearson’s correlation coefficient is not affected by the limitations investigated by Habib et al. (2001) in highly skewed variables. In the cases of the season and distance zones, the sample values of the

_{r}*e*(

*R*) component are determined based on (7) using the overall bias and the deterministic distortion function estimated for that case.

_{r}Figure 9 presents the estimates of the temporal correlation function of *e*(*R _{r}*) at the hourly averaging time scale for the time lags from 15 min up to 3 h. It shows that the errors are correlated in time and that the correlation drops quickly up to the time lags of 60 min. Afterward, it continues to decrease, but more slowly. The quick drop for the time lags below 1 h can be attributed to the nature of the DPA product itself. In fact, as 1-h-running accumulations, the DPA maps separated by time lags below 1 h are not independent by construction because, up to this time lag, the averaging intervals overlap partially. Comparing the three seasons, it is evident that the correlation functions in the random error component for the cold season are significantly higher than those for the warm and hot seasons. On the other hand, the correlations for the warm and hot seasons look similar to each other. However, there is no regular pattern in the dependence of the temporal correlation functions in

*e*(

*R*) on the distance zones. This suggests that, considering the limitations of the information currently available on this question, one curve can represent the obtained temporal correlation functions for all distance zones.

_{r}The estimates of the correlation coefficients of *e*(*R _{r}*) for all the pairs of rain gauge stations at the hourly time scale are presented in Fig. 10 as a function of the distance between the stations. We find that the errors are correlated in space and that the correlation distance is rather large, that is, extending much beyond the RR spatial resolution scale. Also, the scatter in the cold season is smaller and the correlations significantly higher than for the other two seasons. This can probably be explained by the fact that the cold season is characterized by stratiform events that are relatively uniform in space, which is also reflected in the spatial structure of the RR uncertainties. Finally, the presence of the ARS Micronet in zone II is evident from the large number of data points. Unfortunately, this zone is the only one that, thanks to the high density of the Micronet, can provide us with information on the spatial correlations in

*e*(

*R*) at the distances down to a few kilometers. The results for the other zones in Fig. 10 are based on a much smaller number of data points and are probably not very meaningful. However, similar to the temporal correlation functions, no regularity in the dependence on the distance zone is clearly visible for the spatial correlations in the random error component.

_{r}## 6. Discussion

We need to stress that the investigation presented in this study is focused on the exploratory data analysis using mainly the nonparametric functional estimation methods. Although this work was quite extensive and time consuming, it has not exhausted all the issues involved in the PED approach to the RR uncertainty quantification. Among the questions not fully addressed in this study, the most fundamental concern the validation of the presented error model and the obtained estimates, and their transferability to different conditions. The reason that such important problems are beyond the scope of this study is simply their complexity. Their adequate treatment requires additional extensive analyses that we have not performed yet. Below, we briefly discuss the major issues that we plan to investigate thoroughly in our future research projects.

In the most general sense, the validation of a solution (method, model, algorithm, estimate, etc.) is a proof that the solution is compliant with a set of standards that are established to assure its correctness. Except for logically simple and unambiguous problems, there is no universal prescription for these standards that could be used in every situation. In many real data problems, even the seemingly simple notion of “correctness” is not dichotomic and has little meaning without specification of the required accuracy levels. In our case, the validation can also concern different aspects of the results like their sampling error bounds, or their degree of transferability to different situations discussed below. In addition, for a complex multistage analysis such as the one presented here, its validation can concern its different parts from the methodology, through the estimation schemes, to the final numerical outcomes. Therefore, a complete validation of all the procedures and outcomes in this study is a difficult task that can only be accomplished within a new extensive research project.

The first basic validation question in this study pertains to the validity of the statistical-functional representation defined by Eq. (4). The strict definitions of the model components and the straightforward algorithms for their estimation described in section 2 indicate that this representation is a valid and tractable mathematical construct. However, the rigorous formal proof of this assumption has not been obtained yet.

The next important validation question concerns the statistical robustness of the estimates (both parametric and nonparametric) obtained in this study. Although the size of our sample is relatively large, we cannot claim that the sampling errors in the estimated model components are negligible. After thorough consideration of the problem, we concluded that the best way to obtain a reliable quantitative insight into this issue is through computing the error bounds of the estimates. It can be achieved by using resampling methods, such as the bootstrap, to stay within the nonparametric frameworks of this study. The resampling could also reveal possible collinearities between the parameters in the analytical approximations (13) and (15) of our RR error model components. However, this would be a very computationally demanding task, and currently we do not have the resources necessary to perform it within the scope of this study. Some measure of confidence in the hypothesis that the sampling errors in our estimates are not large can be provided by the fairly regular behavior of these results with the changing time scales, distances, and seasons shown in section 4. Nevertheless, this statistical hypothesis needs to be systematically investigated.

A series of validation questions related to each other pertain to the generality of our model and the obtained estimates, as well as to the degree of their transferability to numerous different scenarios not considered in this study. These various conditions have to include (a) different parameter values of the PPS algorithm, (b) other spatial resolution scales and the averaging intervals shorter than 1 h, (c) other radar systems and RR estimation algorithms, and (d) different geographic regions and climatic regimes. Obviously, we are aware that the specific numbers estimated in this study can depend on these situations, just as they depend on the seasons, time scales, and distances considered here. However, based on the experience of this study, we can pose a hypothesis that the general features and qualitative behavior of the error model components will hold at least in most of the new scenarios enumerated above. Another of our working hypotheses is that the dependence of the specific estimates on different conditions will have a regular character that we will be able to interpret and quantify. Of course, speculating about the possible differences between Oklahoma and Florida and Japan, for example, would be fruitless at present. Extensive empirical investigation effort, based on large data samples covering the broad range of the above-mentioned scenarios, is necessary for these hypotheses to be adequately tested.

Some of the questions posed above can be investigated based on the same data sample that we used in this study. The simplest, albeit time consuming, extension of this study elaborates on the effects of different parameter values and correction setups in the PPS algorithm. We will report these results in our future communications. In addition, in our ongoing research, we are working on extending our results from the hourly resolution of the operational WSR-88D products to the shorter averaging intervals of 5 and 15 min. However, at these short time scales, we have to account for the inadequacy of the GR based on single rain gauges. As shown by Ciach and Krajewski (2006), the spatial correlations in rainfall decrease rapidly with decreasing time scales. Therefore, for such short intervals, the point-area differences cannot be neglected in the estimation of RR uncertainties. In addition, as shown in Ciach (2003), at the 5-min resolution the rain gauge data that are processed with the typical tip-counting scheme are significantly affected by tipping-bucket sampling errors. Thus, the spatial rainfall variability becomes an acute problem in RR uncertainty estimation for the shorter averaging intervals. The same difficulty applies also to larger spatial averaging domains. However, efficient methods to account for the GR errors at different spatiotemporal scales based on the available data need still to be developed.

Finally, one of the problems that could not be explored in the present study in all its aspects is the dependence of the estimated uncertainty characteristics on the rainfall regime. Stratifying the dataset according to the three seasons that we used here can be viewed as a poor substitute for a more nuanced approach. It is possible that some appropriate characterization of the regimes on the level of single events, or even the distinctly different parts of the events, could result in the narrower RR error bounds associated with more specific situations. However, we still do not have efficient methods for such a regime-oriented classification, especially for long-lasting and complex events. Therefore, at this stage, we used the simplified seasonal partitioning based on the prevailing surface temperatures. More detailed classification schemes based on different observables (such as current synoptic conditions, surface temperatures, radar echo morphology, etc.), and their possible effects on the RR uncertainty components have still to be investigated.

In our future projects we plan to continue all the extensive analyses discussed above, provided adequate support for these research efforts is available to us.

## 7. Conclusions

We presented a radar-rainfall (RR) uncertainty model that can be applied to the probabilistic quantitative precipitation estimation (PQPE) based on weather radar data. The model has the form of a functional-statistical relationship between RR and the corresponding true rainfall and consists of two components: a deterministic distortion function and a random error factor. The components are estimated based on a 6-yr-long sample of standard WSR-88D rainfall products and the corresponding rain gauge data. The analyses are performed for three seasons, five distance zones, and four time scales. The main findings in this study can be briefly enumerated as follows:

The estimates of overall bias show that the standard WSR-88D rainfall products tend to overestimate the total rainfall in the warm and hot seasons and to underestimate it at distances above 100 km from the radar in the cold season.

After removing the overall bias, considerable conditional bias dependent on RR values still remains, indicating the necessity to account for this systematic effect in a complete RR uncertainty model.

The standard deviation of the stochastic component in the uncertainty model defined by Eq. (4) is typically a monotonically decreasing function of RR value and converges to an asymptote for the largest values of RR.

The frequency distributions of the stochastic component can be approximated with fairly good accuracy by the Gaussian model (in applications, it has to be truncated at zero to avoid negative values).

The stochastic component is correlated in space and time, and the levels of the estimated correlation functions are much higher in the cold season than in the warm and hot seasons.

Comprehensive quantitative evaluation of the statistical robustness, generality, and transferability of our results requires a number of extensive analyses that should be continued in our future research projects.

The results of this study show that the components of the proposed functional-statistical error representation can be fairly well approximated using simple analytical models. The deterministic distortion function can be parameterized by a two-parameter power-law function, whereas the standard deviation of the random component can be approximated with a three-parameter hyperbolic function. The loss of accuracy is unavoidable in such approximations, but it is rewarded by the simplicity of the resulting RR uncertainty description. The fact that, in our specific error representation, the random component estimated for the Oklahoma data sample is close to Gaussian-distributed can facilitate its statistical tractability in further analyses. However, the validity of the particular analytical approximations considered in this study needs still to be investigated in a number of different scenarios discussed in the previous section.

Concerning the spatiotemporal correlations in the random component, the obtained results seem too irregular to venture any reliable parameterization at this stage. The quantitative knowledge of spatiotemporal dependences in the RR uncertainties is necessary in the error propagation studies, or in the PQPE applications, for example. It will be explored in more detail in our future research.

The practical utility of the presented model within the emerging PQPE framework is multifaceted. It allows the construction of maps detailing the probability that, given an observed RR value, the true rainfall exceeds an arbitrary threshold over an arbitrary domain within the radar umbrella. Such maps could be very useful to the flash-flood forecasting problem (Georgakakos 2006). Another application is in the generation of the ensembles of the probable true rainfall arrays conditioned on given RR products. They can then be used as the input for the hydrologic ensemble forecasting models. Such a generator is the necessary element of the PQPE methodology that we are developing, and it will be presented in a separate communication. Finally, the model can lead to a new paradigm of using ground-based rainfall estimates with a known error structure in validation of space-based precipitation products. Discussion of these applications is outside of the scope of this work, and we will present them in the studies resulting from our future research projects.

## Acknowledgments

This work was supported by the NSF Grant EAR-0309644 and the contract under the PQPE initiative by the Office of Hydrologic Development of the National Weather Service. We acknowledge many useful discussions with Richard Fulton, David Kitzmiller, and members of the PQPE Advisory Team. The opinions expressed in this work are those of the authors and do not necessarily reflect those of NSF or NOAA or their subagencies.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Grzegorz Ciach, IIHR-Hydroscience & Engineering, The University of Iowa, Iowa City, IA 52242. Email: gciach@engineering.uiowa.edu