Spatial patterns in thunderstorm rainfall events and their coupling with watershed hydrological response

Weather radar systems provide detailed information on spatial rainfall patterns known to play a signiﬁcant role in runoﬀ generation processes. In the current study, we present an innovative approach to exploit spatial rainfall information of air mass thunderstorms and link it with a watershed hydrological model. Observed radar data are decomposed into sets of rain cells conceptualized as circular Gaussian elements and the associated rain cell parameters, namely, location, maximal intensity and decay factor, are input into a hydrological model. Rain cells were retrieved from radar data for several thunderstorms over southern Arizona. Spatial characteristics of the resulting rain ﬁelds were evaluated using data from a dense rain gauge network. For an extreme case study in a semi-arid watershed, rain cells were derived and fed as input into a hydrological model to compute runoﬀ response. A major factor in this event was found to be a single intense rain cell (out of the ﬁve cells decomposed from the storm). The path of this cell near watershed tributaries and toward the outlet enhanced generation of high ﬂow. Furthermore, sensitivity analysis to cell characteristics indicated that peak discharge could be a factor of two higher if the cell was initiated just a few kilometers aside.


Introduction
Complex interactions exist between spatio-temporal structure of rain systems and watershed hydrological response. While this is a long-standing research issue in hydrology (e.g., [42,61,11,46,29,35,52]), a comprehensive study of these interactions requires detailed rainfall data in space and time that were effectively unavailable until recently.
For many years rain gauge networks were the primary source of storm data. However, these networks are typically sparse and the recorded rainfall data do not adequately represent the spatial variability of the storm (e.g., [32]). In recent years, a new source of rainfall data from weather radar systems has become available. This new technology permits a detailed view of the rainstorm over the watershed with high spatial and temporal resolution not previously available. Access to reliable information on different storm characteristics, 0309 such as the location of the storm over the watershed and its structure and direction of motion, presents a new stage in rainfall-runoff analysis.
It was hoped that the use of detailed radar rainfall information as input into hydrological models, representing the watershed hydrological response, would significantly improve understanding of rainfall-runoff processes and aid in predicting their outcome. However, a review of studies using radar rainfall data in hydrological modeling (e.g., [39,24,34,44,2]) provide mix evidence for such improvements [4,49]. Among the difficulties reported were the accuracy of radar rainfall estimates [1,27], and runoff sensitivity to subpixel rainfall variability [25,32,60,22]. Arguably, even though information about rainfall patterns exists in radar rainfall data, there is a need to develop new ways to exploit this information and gain greater insights into rainfall and subsequent watershed response behavior [19,56].
The current study suggests an innovative approach to the above challenge. Our focus is on semi-arid watersheds and their response to air mass thunderstorm events. We emphasize, as a key issue, the spatial structure at which the radar rainfall data are represented in the hydrological models. In the standard approach, radar-based rainfall data are applied into the model in a grid structure that states the amount of rain at each pixel in a given time step. The grid is either that of the observing radar system, or a transformation is performed into another standardized grid. For example, radars in the US National Weather Service (NWS) weather radar network observe rainfall in polar grids with 1°· 1 km pixels for different tilts with a time sampling interval of 5 min. Data from several radars are merged and integrated to produce surface hourly rainfall at national grids of about 4 · 4 km 2 pixels. This radar rainfall product is sent to the River Forecast Centers (RFC) and serves as input to NWS hydrological models [15].
Although grid structures are common and straightforward, some difficulties regarding their use in hydrological modeling are raised here: 1. Often, hydrological modeling is interested not only in simple runoff estimation but also in identifying the key elements affecting runoff production. In such cases, the models are structured to allow questions about interactions between hydrologically significant elements. For example, the surface runoff processes in the watershed are commonly conceptualized as occurring over hillslopes and channels. In these cases, the topographical data are not used in the model in their original structure (grid or contours) but rather, the data are processed to identify the hillslope and channel elements and their characteristics (slope, length, area, etc.). The same argument holds for rainfall data. In many cases, only limited information can be attained by representing the rain as an arbitrary data structure (i.e., a grid of pixels with uniform rainfall). On the other hand, identifying Ôsignificant elementsÕ within the rainfall data, may allow advanced analysis by enhancing our ability to concentrate on the interaction of the rain system with the hydrological units and processes.
2. The representation of rainfall in a grid structure implies a strong assumption concerning uniform rainfall over a pixel. This assumption has important consequences in hydrological modeling, especially in a semiarid environment where runoff generation processes are sensitive to small-scale rainfall variability (e.g., [11,17,18,5]).
3. Radar rainfall data are known to contain relatively large errors resulting from the indirect measurement [27]. To reduce error, radar rainfall data are often smoothed out, usually by integration to larger pixels and time steps. In other words, the grid structure is maintained but with lower resolution. A tradeoff exists between reducing error and losing important information contained in the small-scale data [16].
In the current study we introduced rainfall input data to the hydrological model in the form of hydrologically significant rain elements instead of the standard radar-grid structure. Because the analysis is focused on air mass thunderstorm rain events, the hydrologically significant elements chosen were the basic convective units, usually referred as ''rain cells''. Rain cells are conceptualized in our work as isotropic circular elements with maximal rain intensity at the cell center. The location, maximal intensity and decay factor of each rain cell in the investigated storm were fed into the hydrological model and the watershed response was computed. In this approach, the computed runoff is directly linked to spatial characteristics of the rainfall and thus, we believe, greater insights can be gained on interactions of rain and hydrological systems (see also [36]).
The objectives of this study are 1. To represent rain fields of air mass thunderstorm events as a set of rain cells (the rain cell parameters are derived from observed weather radar data by applying a conceptual rain cell model). 2. To compare the spatial characteristics of the rain fields, as obtained by applying the rain cell model, to the ones obtained by using grided radar or dense gauge network data. 3. To identify the dominant factors of the rain cells by applying the cell-based rain fields of an extreme thunderstorm event to the hydrological model. 4. To examine runoff sensitivity to rain cell characteristics identified in objective 3.

Study areas and data
The two study areas are located in the semi-arid climate regime of southern Arizona. The summer weather of this region is strongly affected by the North American monsoon [8,54], which results in frequent air mass thunderstorms that are highly convective, intense, localized, and of short duration.
The 148 km 2 Walnut Gulch Experimental Watershed (WGEW; [18]) is located 50-70 km east-southeast of the Tucson WSR-88D weather radar, which is a part of the NWS weather radar network (Fig. 1a). The watershed is equipped with a dense network of rainfall gauges (Fig. 1b) managed by the Agriculture Research Service of the US Department of Agriculture (USDA-ARS). Thirteen convective storms from the 1999 and 2000 monsoon seasons (June-August) were selected for the study (Table 1). Radar data from a 1125 km 2 segment encompassing the watershed were used for the analysis (see Fig. 1b). Because radar beams at the first and second tilts are partially blocked by terrain before they reach the watershed, only the third tilt data (elevation The radar segment (azimuth: 93-122°and range: 42-78 km relative to the Tucson radar) encompassing the 148 km 2 WGEW study area and the 74 gauges in the watershed. The radar segment of the Phoenix study area is in azimuth 291-320°and range 42-78 km relative to the Phoenix radar. Radar polar grid resolution is of 1°· 1 km. (c) The radar segment of the Phoenix study area (azimuth 291-320°and range 42-78 km relative to the Phoenix radar) and the gauge network of the Flood Control District of Maricopa County. Radar polar grid resolution is of 1°· 1 km. angle of 2.4°) were used for the analysis (equivalent to 3 km altitude above ground over the study area; average watershed elevation is roughly 1340 m above mean sea level). The relatively high altitude might impose an additional error in the radar rainfall estimates resulting from evaporation or hail contamination. However, it should be emphasized that a very limited portion of the US has NWS radar data available within 2 km of the surface and in the western US radar data are typically sampled near and above 3 km [28]. Rain data from 74 rain gauges located in the watershed area ( Fig. 1b) were used for radar calibration and in the evaluation procedure.
The second study area is the Phoenix metropolitan area, located northwest of the Phoenix NWS WSR-88D weather radar (Fig. 1a). Data from a summer storm in 2002 were used for this study ( Table 1). The radar data segment used to cover the Phoenix storm was similar in size and distance (from radar) to that used for the WGEW events. Because there are no topographic blockage problems in the Phoenix area, we were able to use radar data from the lowest three tilts (roughly equivalent to 1, 2, and 3 km above ground level) for the analysis. Gauge data for the Phoenix area were used for examination of the derived radar rainfall estimation equation (see Section 2.2 below). Thirty-four gauges operated by the Flood Control District of Maricopa County are within the analyzed segment of the radar (Fig. 1c).

Radar rainfall estimation
Radar reflectivity data (Z) [mm 6 m À3 ] are represented in polar coordinates centered at the radar station at a resolution of 1°· 1 km (equivalent to 0.88-1.23 km 2 at the examined areas) for different tilts with a time sampling interval of 5 min. Reflectivity data are commonly converted to rain intensity data (R) [mm/h] using a power law Z-R relationship Z = aR b . In this study, the exponent parameter (b) was set to the value of 1.4 that is used by the NWS for convective rainfall [15], while the multiplicative parameter (a) was adjusted based on comparison of gauge and radar storm depth data. The resulting Z-R relationship, based on analysis of radar and gauge data for the selected 13 storms over WGEW, was An upper threshold of 100 mm/h was applied to the estimated rain intensity. This is a default threshold rain intensity used by the NWS to reduce unreasonably large estimates caused by hail cores in thunderstorms [15]. The Z threshold value used here (56 dbz, where 1 dbz = 10 Log Z) is different from the threshold used by the NWS (53 dbz, [15]), because of the different Z-R relationships.
The above Z-R was examined for the Phoenix storm by comparing total storm depth of 34 gauges within the analyzed segment to radar estimates in pixels above the gauges. The total bias found in radar estimates was À9%, +3% and À22% for the first, second and third tilt data, respectively. Further details pertaining to radar rainfall estimation can be found in Morin et al. [38].
3. Representing spatial rainfall patterns in the form of rain cells

Rain cell model
Several mathematical descriptions were suggested for rain cell shapes. Circular uniform cell shape was assumed in the stochastic modeling studies of Eagleson [9], Cox and Isham [6], Wheater et al. [58] and Marco and Valdes [29]. Uniform, exponential decay and quadratic exponential decay functions for circular rain cells were examined in the stochastic modeling studies of Rodriguez-Iturbe et al. [45] and Eagleson et al. [10]. The latter found the quadratic exponential function most suitable for representing spatial variability of storm rain depth over the WGEW using rain gauge. More complex descriptions of rain cell shapes use uniform, exponential decay function, quadratic exponential function or their combinations, in a bivariate form to get an elliptic cell shape [3,40,13,59,53]. In most of the above studies, mathematical representation of rain cells is part of stochastic rainfall modeling with the motivation to randomly generate simulated rain fields that are similar to observed rain fields in the statistical sense. In the current study we use the mathematical cell shape description to decompose observed radar rain data into set of rain cells and replace the observed radargrid data with rain fields induced by the retrieved rain cells. For this goal we adopted the circular quadratic exponential shape that was already tested for the studied area [10].
The rain fields are assumed to be composed of multiple rain cells with the rain intensity field at rain cell i is given by where R i (d) is rain intensity [mm/h] for rain cell i at a distance d [km] from the rain cell center at coordinates (X i , Y i ), b i is the rain intensity [mm/h] at the center of rain cell i, and, a i is the decay parameter [km À1 ] of rain cell i. Eq. (2) defines a two-dimensional non-normalized Gaussian surface, where b i represents the amplitude and a i represents the spatial extent. Parameter a i is equivalent to the inverse of two times the standard deviation of the Gaussian distribution. Accordingly, relationships between the rain cell model parameters and other rain cell characteristics can be derived: 1. If I i is the integral of rain intensity over rain cell i [km 2 mm/h], then Rearrangement of Eq. (4) and applying it to two different threshold values s 1 and s 2 results in Eqs. (3) and (5) are used in the rain cell retrieval algorithm, explained in the next section, to estimate parameters b i and a i , where I i and A s i are derived from the observed data.
The selected rain cell model maintains a relatively simple structure while still representing hydrologically important spatial characteristics such as: location and magnitude of maximum rainfall intensity, rain-no rain areas (rain cell coverage) and small-scale variability (within rain cells). As aforesaid, the model validity for the studied region is supported by the work of Eagleson et al. [10]. It should be emphasized however, that the circular shape is not expected to be applicable for rain systems containing un-isotropic elements, such as fronts or other linear shape systems. Whether this rain cell model is valid for air mass thunderstorms in other locations still needs to be examined.

Rain cell retrieval algorithm
The following is a brief description of the procedure used for retrieving rain cells for each radar map (the entire algorithm is available from the first author on request). A radar map represents a rainfall value (reflectivity or rain intensity) for each polar pixel in the analyzed region at a given time step. An example illustrating the computational stages is presented in Fig. 2.
Step 1. Input: radar reflectivity map for a specific time step (Fig. 2a).
Step 2. Derive radar rain intensity map by converting radar reflectivity data at each pixel into intensity using the Z-R relationship described in Section 2.2 ( Fig. 2b).
Step 3. Divide the rain map into individual segments, where each segment is defined as the region dominated by a single rain cell (Fig. 2c). A segment is identified by a local maximum rainfall pixel that is not part of any existing segment. The segments are gradually expanded to neighboring pixels with rainfall values lower than or equal to the ones already existing in the segment but higher than a rainfall threshold (5 mm/h, user defined parameter). The search for local maximum pixels is from high to low values. This ensures that pixels with more than one segment affiliation are assigned to the segment having higher maximum value.
Step 4. Eliminate or merge segments with small area or low rain intensity. After identifying all segments in a rainfall map, segments with area smaller than a minimum threshold (9 km 2 , user defined parameter) or maximum pixel rainfall lower than a threshold center pixel intensity (25 mm/h, user defined parameter) are either eliminated (if isolated, i.e., no other segment is neighboring) or merged with the nearest neighbor segment.
Step 5. Identify rain cell for each segment by extracting rain cell parameters. For each segment the center location of the maximum pixel is considered the rain cell center (black points in Fig. 2c). For each rain cell i, the integral of rain intensities over the rain cell, I i , is estimated as the sum of rain intensities over the rain cell data segment. Rain areas greater than the threshold, A s i , are estimated as the total area of the data segment pixels with rain intensity higher than a given threshold for a series of rain intensities (in the range of 5-100 mm/h). The model parameters are estimated using I i and A s i . a i is estimated using Eq. (5) (Fig. 2d) and b i is then calculated using Eq. (3). The estimated b i value can, in principle, be lower than the threshold cell center intensity (step 4) or higher than the hail rain intensity upper threshold (see Section 2.2).
Step 6. Output: number of rain cells for each radar map and the four model parameters for each cell (Fig. 2e).
Step 7. Perform steps 1-6 for the series of radar maps of a given storm through time. The algorithm produces a description of the rain cells (number, location, and parameters) as they evolve throughout the storm. The threshold parameters used in the algorithm discriminate between pixels that comply with the assumed rain cell structure and pixels that do not. The threshold values selected for the current analysis (i.e., area of rainfall above 5 mm/h is at least 9 km 2 and rain intensity at the maximum pixel is at least 25 mm/h) are within the range suggested by other authors (see for example [52] and references therein). Additionally, results are also presented for the case where the maximum pixel intensity threshold is lowered to 10 mm/h.
The above algorithm assumes no significant interacting area between adjacent rain cells. Therefore, we can segment the data and estimate parameters for each rain cell independently. This assumption appears reasonable for the current air mass thunderstorm dataset. However, for a more general case, a global optimization procedure may be necessary in order to estimate the parameters of multiple rain cells simultaneously.

Retrieved rain cell characteristics
The algorithm was applied to radar data from 13 storms over the WGEW (Fig. 1b, Table 1). Statistics of rain cell characteristics is presented in Fig. 3 and Table 2. As an example, we examine the results for the July 6, 1999 storm over WGEW (Fig. 4). The intense phase of the storm lasted roughly 1.5 h (time steps [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Within this period the number of rain cells in the study area at each time step ranges from 3 to 10 and the b and a parameters range from 10 to 200 mm/h and from 0.2 to 0.7 km À1 , respectively. After this high intensity period, the number of rain cells was reduced to 0-5 cells, with b ranging from 10 to 30 mm/h and a ranging from 0.15 to 0.45 km À1 . If the threshold for center pixel intensity (see Section 3.1, step 4) is reduced from 25 to 10 mm/h, the second storm period (time step 20 and on) is characterized by a large number of low intensity rain cells (Fig. 4b-d). These differences at the dissipating stage of the storm suggest that the convective rain cells are no longer the governing structure and widespread rainfall is more dominant in the rain field. Under these conditions the assumption of circular rain cell elements does not necessarily hold and the widespread rainfall is interpreted as many small and low intensity rain cells. Keeping the threshold for center pixel intensity at its high value (25 mm/h) prevents detection of these artificial rain cells.
The algorithm was also applied to radar data for one storm over the Phoenix region using the same Z-R relationship and threshold parameters used for the WGEW study area. Rain cell characteristics at different radar tilts were compared for the Phoenix data ( Table 2, Fig. 5). The first tilt data generated more rain cells, including times at which no rain cells were identified at the second and third tilts. Close examination of the rain maps reveals that in most cases these cells contain ground clutter contamination at the first tilt, while the others appear at the dissipation stage of the storm. Comparison of the second and third tilt data showed a relatively high level of matching rain cells (in terms of their number and parameters). At time steps 15-18 for the second tilt there is a considerably larger number of rain cells relative to the third tilt. As suggested above, the difference in rain cell number can probably be attributed to the dissipating stage of the storm occurring at these time steps. In addition, at time step 10, a rain cell with b = 218 mm/h was found at the third tilt while the associated rain cells at the second and first tilts had considerably lower b values of 148 mm/h and 104 mm/h, respectively. However, these high reflectivity data (high-er than 63 dbz), probably indicate hail contamination, in all three tilts. In the third tilt, the high reflectivities covered a larger area, which resulted in a higher integral   )). Except for this case, the range of the b and a parameters was similar for the second and third tilts. Although further analysis is needed for fine characterization of the rain cells along the vertical profile the analysis conducted over the Phoenix area suggests that the third tilt data are both suitable and sufficient for the analysis presented here. We emphasize this point because, as already mentioned, the analysis of the WGEW storms used only the third tilt data because the first two tilts are blocked by terrain.

Evaluating spatial characteristics of cell-based rain fields
An important step in the study was to evaluate the spatial characteristics of the rain fields as obtained by the rain cell model. This step includes comparisons of these spatial characteristics with those of the grided radar data and of rain fields computed from the dense gauge network in the WGEW (an average of one gauge in 2 km 2 ). Although the gauge data had already been used to adjust the Z-R relationship and thus to remove the overall bias from the radar rainfall estimations (see Section 2.2) they still provided a set of relatively independent observations. Gauge rainfall maps, with a 100 m pixel size, were generated from 3-min gauge-rain intensities, using the multiquadric interpolation method [52]. A time delay of 5 min was applied to gauge data to account for the altitude sampling differences [37]. A rain gauge derived rainfall map was generated at each time step for which radar rainfall maps were available. Since gauge rainfall observations are available only within the watershed area, radar-based rain intensities outside of the watershed boundaries were assigned to zero.
The evaluation process includes comparison of several rainfall sources as follows: 1. Rain cells derived from radar data, hereafter called radar-rain cells. 2. Rain cells derived from rain gauge data are termed gauge-rain cells. The interpolated gauge rainfall maps were used as input to the retrieval algorithm described above (Section 3.1 start at step 3). In order to use the exact same procedure and avoid scale related discrepancies, the rain fields were remapped first to the radar polar grid. Also, the same threshold parameters were used in the algorithm for the two inputs. 3. Observed radar rainfall data at the radar-grid, hereafter called radar-grid. 4. Interpolated gauge rainfall, hereafter called gaugeinterpolated.
Below we describe several evaluations tests that were applied in order to examine the different rain field spatial characteristics (also see Fig. 6).

Direct comparison
In this test, the gauge-rain cells and radar-rain cells over the WGEW were compared. Gauge-and radar-rain cells were matched according to the distance between the cells, allowing up to 4-km difference. Seventy-one percent of the gauge-rain cells were identified from the radar data. In contrast, 40% of the radar-rain cells did not have a matching cell in the gauge data. The correlation between the derived parameters for 398 matched rain cells is r = 0.65 for the b parameter and r = 0.50 for the a parameter. Over the matched rain cells with distances up to one kilometer the correlation of the b and a parameters is 0.81 and 0.51, respectively. Estimations of both parameters are higher for radar-rain cells (27% and 21% for b and a, respectively) due to the following possible reasons: (1) Variations in rain intensity from rain aloft to ground rain (due to evaporation for example) can cause higher b values for radar-rain cells relative to gauge-rain cells. It may also cause differences in cells number if rain exceeds the intensity threshold aloft but does not on the ground.
(2) Although bias between accumulated rain depth in gauges and in radar pixels above gauges was removed in the calibration process (Section 2.2), it still possible to have bias between gauge and radar rainfall over the entire watershed. As shown later (Section 4.4) a positive bias is found for the radar data, resulting in higher values of the b parameter for radar-rain cells.
(3) The average gauge density in WGEW is one gauge per two km 2 (roughly two radar pixels). This causes gauge-rain field to be smoother than the radar field, which results in lower estimations of a.
(4) A certain number of rain cells were located close to the watershed divide. When retrieved from gauge data only partial rainfall information is available, causing them to be either missed or detected with large error in parameters.
(5) Timing discrepancies due to the different altitude (3 km) of radar and gauge data depend on the dropsÕ terminal velocity. For the same group of storms, Morin et al. [37] found a time delay of 3-9 min between radar and gauge observations. A 5-min time delay was applied to the gauge data to roughly represent averaged conditions. However, the actual time delay varies in space and time (see for example [20]) and can result in a mismatch between the gauge and radar data.
(6) Space discrepancies depend on the vertical profile of the horizontal wind velocity and can lead to several kilometers distance between radar-derived and gaugederived rainfall location. We partially overcame this problem by allowing up to 4 km distance (see [37]) between matched gauge and radar derived rain cells in the direct comparison test. However, since the actual space shift varies, spatial discrepancies can still affect the direct comparison evaluation results.
In summary, because the direct comparison test examines the retrieved rain cell parameters, it is relatively sensitive to gauge-radar data discrepancies. The following three evaluation tests do not match individual cells, and are less sensitive to these types of errors.

Parameter distributions comparison
Here, the gauge-and radar-rain cells are compared in terms of their parameter distributions. Fig. 7 presents the parameter distributions for radar and gauge-rain cells within the watershed. In general, the distributions are similar in their shape but somewhat different in their central tendencies. The b and a distributions of radarrain cells trend toward higher values. In addition, more radar-rain cells have b values lower than 30 mm/h. Some of these rain cells probably did not pass the threshold level on the ground (only 30% of cells had a matched gauge-rain cell). The distribution of cell areas (area above threshold rain intensity of 5 mm/h) indicates larger gauge-rain cells (Fig. 7d), suggesting a larger sampling distance of the gauge data (see Section 4.1).

Spatial correlation comparison
Spatial correlation is the correlation of rainfall at two points separated by a given distance. Spatial correlation curves describe the decay of correlation with increasing separation distance and are often used to characterize the spatial variability of rainfall patterns [63]. Here, we compared the spatial correlation curves of the following rainfall fields: gauge-interpolated, radar-grid, gauge-rain cells, and radar-rain cells. In addition, we calculated correlation curves based on radar-rain cells with a 20% difference in their a parameter. In the analysis, the correlations were calculated based on rain intensity at points of gauge locations (in 0.5 km separation distance bins). We included all time steps over the 13 storms in Fig. 6. Evaluation scheme. The four rainfall fields used in the evaluation process and the four evaluation tests. Arrows indicate the rainfall fields used in each test. which one or more gauges recorded a rain intensity higher than 25 mm/h. The comparison of the spatial correlation curves is shown in Fig. 8.
The rain field induced by gauge-rain cells exhibits higher spatial correlations as compared to that of radar-rain cells and fits the À20% a curve. The differences in spatial correlations are probably due to the differences in spatial sampling between gauge and radar that result in smoother gauge-based rain fields. At short separation distances, spatial correlations of rain fields induced by gauge-rain cells and radar-rain cells are higher than those of gauge-interpolated and radar-grid rain fields. This is explained by the effect of within rain cell organization that is dominant at short distances and results in smoother rain fields. At larger distances, the correlation of the rain cell field drops rapidly and becomes lower than the non-rain cells correlations.

Areal rainfall comparison
This test estimated the correlation and bias of rain intensities between gauge-rain cells, radar-rain cells and radar-grid to gauge-interpolated rainfall. This was done by averaging rain intensity over the watershed and in sub-regions of 4 · 4 km 2 and 1 · 1 km 2 (total of 11 averaging areas; Fig. 9a). As expected, the highest correlation with gauge-interpolated rainfall was found for gauge-rain cells. At the watershed scale, very small and insignificant difference exists between correlations of radar-grid rainfall with gauge-interpolated rainfall and of radar-rain cells with the gauge-interpolated rainfall. Over the sub-regions the latter is significantly smaller in the range of 4-16%. Rainfall of gauge-rain cells and radar-rain cells have large negative biases resulting from applying thresholds for rain cell pixels (5 mm/h) and cell center intensity (25 mm/h). Reduction of the center pixel intensity (see Section 3.1, step 4) from 25 to 10 mm/h results in smaller negative bias (Fig. 9d). Fig. 9 also indicates a positive bias of the radar-grid rainfall relative to the gauge-interpolated rainfall. This bias includes radar pixels without gauges as opposed   8. Spatial correlation curves of the following rainfall fields: gaugeinterpolated, radar-grid, gauge-rain cells, radar-rain cells, radar-rain cells with À20% difference in their a parameter and radar-rain cells with +20% difference in their a parameter. The correlations were calculated based on rain intensity at the points of gauge locations in 0.5 km separation distance bins, including all time steps over the 13 storms at which one or more gauges recorded rain intensity larger than 25 mm/h. to the bias computed in the calibration process (Section 2.2).

Summary of evaluation results
1. Application of the rain cell retrieval algorithm to radar data identified most (71%) of rain cells identified from gauge data at a relatively close distance (2 km on average). A moderate correlation between radar-rain cells and gauge-rain cells parameters was found (0.65 and 0.5 for the b and a parameters, respectively).
2. Radar-rain cells are characterized by higher (+27% bias) b parameter values relative to gauge-rain cells.
Possible explanations include altitude differences, the effects of evaporation and a positive bias in radar rainfall relative to gauge rainfall for the entire watershed.
3. Gauge-rain cells are in general larger (+52%) than the radar-rain cells and the rain fields are characterized by higher spatial correlation (Fig. 8). We suspect this to be the result of inadequate spatial sampling by the gauge network. The outcomes are lower a parameter values (À15%) and smoother rain fields.
4. Rain fields induced by rain cells (from gauge or radar data) represent only rain pixels that pass the selected thresholds. This results in an underestimation (À32% on average) of areal rainfall as compared to the gaugeinterpolated rainfall. Correlation of areal rainfall between radar-rain cells and the gauge-interpolated rainfall is relatively good (0.89, 0.75, 0.66 for the entire watershed, the 4 · 4 km areas and the 1 · 1 km areas, respectively).

Cell-based rain fields as input into hydrological model
The above sections presented retrieval of rain cells from radar data. Each rainfall map is translated into a set of rain cells specified by their location, maximum rain intensity, and decay factor. The rain cells are hydrologically significant elements and may play an important role in the process of runoff generation. In this section we used rainfall in the form of rain cells as input to a hydrological model to examine the watershed hydrological response in an extreme case study.
The hydrological model used is KINEROS2. It is a physically based, event oriented, rainfall-runoff model [62,50] developed by USDA-ARS scientists for watersheds in semi-arid environments. The model represents the watershed as a cascade of overland flow planes and channels, thereby allowing rainfall, infiltration, runoff and erosion parameters to vary spatially. Recently, a GIS-based tool (AGWA-see: www.tucson.ars.ag.gov/ agwa) was developed by the USDA-ARS [33] for delineating watersheds into hillslope contributing areas (abstracted into overland flow plane model elements) and channels. AGWA also generates model parameter files, based on topography, soil and land cover information. In this study, we used the KINEROS2 model with default parameters generated by AGWA for the WGEW (delineation of 53 planes-average area 2.8 km 2 , and 21 channels). Initial soil moisture content was estimated from daily rain depth data using simulation model [21,31]. The study does not include sensitivity analysis to watershed characteristics and therefore all parameters were kept as their default values throughout the analysis.
We analyzed the extreme rainfall-runoff event of August 11, 2000, which totaled to 25 mm watershed Fig. 9. Comparison of areal rain intensity from radar-grid rainfall, gauge-cell rainfall and radar-cell rainfall with gauge-interpolated rainfall. The comparison is based on all 13 storms and at 11 locations and spatial scales: the entire watershed area, five 4 · 4 km areas (L1-L5) and five 1 · 1 km areas (S1-S5). (a) Evaluation sites and gauge locations, (b) correlation, (c) bias, and, (d) bias with 10 mm/h threshold for center pixel intensity (reduced from 25 mm/h). average rainfall depth with a maximum gauge depth of 91 mm. Most of the rainfall occurred in less than an hour and characterized by extreme rain intensities (e.g., maximal 1-min gauge intensity of 391 mm/h). The thunderstorm caused a high runoff flow in the watershed with a peak discharge of 154 m 3 /s at the outlet, the highest recorded peak since 1957.

Storm decomposition into rain cells
Two radar-based rainfall inputs were analyzed using the hydrological model, the grid rainfall data and the rain cell data. The rain cell data were generated by applying the retrieval algorithm above to the storm radar data. Fig. 10 presents comparison of observed vs. computed runoff in terms of peak discharge rates at a range of subwatershed scales (Fig. 10a) and hydrographs at the watershed outlet (Fig. 10b). It should be emphasized that the radar rainfall estimates were based on the Z-R relationship described at Section 2.2. If the default NWS Z-R relationship for convective precipita-tion is used [15], the computed runoff peak discharge is more than two times higher [38]. The computed runoff hydrographs based on the two rainfall inputs are reasonably close to each other and to the observed runoff in terms of peak discharge. In terms of time to peak, the computed runoff is on average 25 min early relative to the observed runoff (not shown). This difference is equivalent to 16% error, which is reasonable considering the complexity of rainfall-runoff modeling in semi-arid watersheds [32].
Comparing the two computed runoff hydrographs presented in Fig. 10 revealed relatively small differences. This suggests that both inputs contain essentially the same information required for predicting the hydrological response (as represented by the hydrological model). However, we believe that the rain cell input also allows the characterization of the major storm factors in runoff generation.
A manual examination and tracking of the rain cells in the August 11, 2000 rainfall-runoff event (Fig. 11, Table 3) identified five rain cells (cells A-E) that developed close to or over the watershed. Two of the rain cells (A and C) were intense (in terms of maximum intensity and volume) and lasted a relatively long time (more than 60 min). Rain cell A initiated outside and north of the watershed and then moved southward into the watershed. Rain cell C developed within the watershed and moved west-northwest toward and beyond the watershed outlet. The three other cells (B, D and E) moved generally northward, had shorter duration (15 min) and were less intense.
Comparison of the computed runoff response for each rain cell separately indicates that only cell C generates a significant peak discharge at the watershed outlet (Table 3). This is probably due to its location within the watershed for most of its life cycle. Another important factor is the specific configuration of rain cell C relative to the watershed. Examining the detailed response of each model element (overland flow planes and channels) revealed that the cell passed close to watershed  tributaries at three locations, precipitating more than 25 mm of intense rainfall over their associated contributing areas. This generated significant excess rainfall over the associated runoff model elements and high peak flows at the tributaries to the main channel (Fig. 12). The flow towards and along the main channel was such that runoff peaks arriving from the tributaries were relatively close in time to the runoff peak of the main channel, resulting in an intensification of the peak flow. Over the last 7 km of the channel, there was no lateral contribution to the flow from hillslope areas or tributaries, which resulted in reduction of peak flow due to channel transmission losses.

Sensitivity analysis
Sensitivity of watershed rain depth and computed peak discharge to maximum cell intensity and its spread was examined by applying a factor (0.5-1.5) to the b and a parameters of cell C (Fig. 13a). Rainfall over the watershed increases linearly with b and decreases exponentially with a, as dictated by Eq. (2). Changes in rainfall are directly and non-linearly translated into changes in runoff peak discharge (Fig. 13a).
Sensitivity to both parameters is relatively high with somewhat higher sensitivity to the a parameter. For example, decreasing b by 10% results in a decreasing of watershed rain depth by 10% and a decreasing of 35% in runoff peak discharge. Decreasing a by 10% increased watershed rain depth and runoff peak discharge by 20% and 65%, respectively. The functional relationships of runoff peak discharge and watershed rainfall are presented in Fig. 13b.
Sensitivity to rain cell location was examined in two ways: cell direction and cell starting point. Sensitivity to cell direction was examined through altering the locations of rain cell C relative to its most southeastern point (Fig. 14). The two peaks seen in the rainfall graph reflects the two major directions (north-west and northeast) at which the rain cell precipitates rainfall over the watershed for the longer period of time (see Fig. 11). Runoff peak discharge is maximized 10-15°c lockwise to direction of maximal rainfall. These preferable directions are dictated mainly by amount of rainfall  12. Computed response of model elements in relation to rain cell location and trajectories for the 8/11/2000 storm over the WGEW. Averaged rain depth from the modeled rain cell C (see Fig. 11) to each overland flow plane is shown. For selected points along the channel network the time of peak discharge (T p ) and peak discharge (Q p ) are presented.
precipitated over the watershed and secondarily by the closeness of the cell track to the main channel. Note that simulated peak discharge at the actual direction of the rain cell (0°change in Fig. 14) is very close to its maximal value for different directions. On the other hand, in terms of starting point, computed runoff peak increases more than two fold with less than a 4 km shift in cell location (Fig. 15). The figure presents watershed rainfall (Fig. 15a) and outlet runoff peak discharge (Fig. 15b) for different starting points of cell C. As opposed to the above sensitivity tests, where watershed rainfall amounts were the most important factors in determining runoff peak discharge (e.g., Fig. 13b), in the current test, cell track position relative to the main channel is a major factor as well.

Discussion
This paper deals with representation of rainfall spatial patterns associated with air mass thunderstorm events in hydrological model. The rainstorm is represented as a set of convective rain cells retrieved from radar data by application of a conceptual rain cell model. The rain cells characteristics are utilized as input to the hydrological model. These characteristics include all main spatial features of rainfall patterns: location and magnitude of maximum rainfall intensity, rain-no rain areas (rain cell coverage) and small-scale variability (within rain cells). These features are important characteristics of the rain system and are known to play a significant role in watershed response to rainfall. We believe that through their direct linking to the hydrological model, new insights can be acquired in understanding the behavior of the rain and hydrological systems and their interaction.
Fourteen air mass thunderstorm events over the Walnut Gulch Experimental Watershed and the Phoenix area in Arizona were analyzed. The rain cells center rain intensity parameter (b) varied over a relatively large range (mean = 61.5 mm/h, standard deviation = 42.3 mm/h), while the distribution of the areal spread parameter (a) showed small variability (mean = 0.39 km À1 , standard deviation = 0.08 km À1 ). Physically, these observations indicate that at 2.5 km (2.1-3.2 km) distance from the cell   center the rain intensity is reduced to 14% of its intensity at the center. Rapid decline in rain intensity was previously reported in other studies of the Arizona region [43]. Our estimation of a is supported by Eagleson et al. [10], which reported an estimated mean value of a = 0.43 km À1 and standard deviation of 0.17 km À1 over a set of 426 storms. The referred work used a similar rain cell model, but in a stochastic framework, for WGEW gauge storm depth data. The objective was to produce similar spatial distributions of total storm depth from the model and the observed data. The relatively stable and narrow range of a suggests its generality, possibly as a climatic characteristic representing the areal extent and spatial variability for air mass thunderstorms. This assumption can be examined by applying the model to radar data from thunderstorm rain events in other regimes.
We explored the potential use of rainfall, represented as rain cells, as input to a distributed hydrological model. Although the original radar rainfall data produced a relatively similar computed runoff hydrograph as the rain cells data, it is the added insights on the spatial watershed response that makes this approach beneficial. The retrieval algorithm process enables the decomposition of complex rainfall patterns into rain cells with defined parameters. The explicit representation of rainfall spatial patterns in hydrological model input allows derivation of a more comprehensive link between runoff response and spatial rainfall patterns. Reviewing numerous studies that investigated the relationship between the two, revealed the general use of two approaches: (1) Synthetic rainfall is used to simulate different spatial patterns (e.g., [51,26,41]); or, (2) Spatial characteristics of observed rainfall data are statistically related to the associated runoff response [30,52]. Using a rain cell representation of ''real storm'' data allows a combination of the two approaches. As in the first approach, a spatial structure is assumed (the rain cell shape), representing our a priori understanding, but the parameters of the assumed structure are derived from observed data and thus, as in the second approach, real observations are integrated into the analysis.
While Syed et al. [52] employed the same model for describing rain cells they only used this model to evaluate several interpolation algorithms. This was accomplished by superimposing a theoretical rain cell produced using Eqs. (3) and (4) over the center of the Walnut Gulch Experimental Watershed; sampling the theoretic rainfall at actual rain gauge locations, using a number of interpolation schemes to produce a 100 · 100 m grid of the rainfall; and, evaluating how well the interpolated rainfall field reproduced the theoretic field. After that, the selected interpolation approach (multiquadric), was only used to interpolate rain gauge storm observations onto a 100 m grid from which storm characteristics were numerically computed (centroid, volume, distance of watershed outlet, etc.). These observed storm ''metrics'' were regressed onto runoff observations to investigate how well they could statistically predict watershed runoff.
For the research described herein it is important to understand the implications of using grided rainfall when comparing the runoff response to the two rainfall inputs (rain cells and grid). Grayson and Bloschl [19] pointed out that a uniform assumption imposes a certain spatial organization to the elementary unit structure. Assuming uniform rain intensity over a pixel (the elementary unit in this case) may considerably affect runoff response, in particular for semi-arid environments where runoff to rainfall ratios are small [11,18,5]. To retain small-scale rainfall variability of rainfall input, the conventional approach uses high-resolution ''noisy'' rainfall data. In contrast, the spatial fields induced by a set of rain cells are smoothed by the procedure of Gaussian surface fitting while still allow representation of sub-pixel variability using the continuous mathematical function. It should be emphasized however, that although rainfall variability can be represented at small scales as we wish, it cannot be confirmed without rain gauge data with spacing at that scale or smaller.
When sufficiently long records of rainfall data are available, the distributions of the rain cell characteristics can be analyzed [12,53]. This will allow examination of the return-period extremity in terms of the rainfall spatial structure as compared to the conventional approach using the frequency of rain accumulation for a given duration (see for example, [7]). These distributions could then be derived for ungauged regions assuming reliable radar rainfall data are available. If corresponding runoff data are also available, the extremity of the rainfall event can be related to the associated runoff response directly or statistically (e.g., [48,64]). Such analyses have the potential to provide new insights in estimating the ''probable maximum precipitation'', ''probable maximum flood'' and other related extremes of interest [14,47,29].
Finally, the focus of this paper has been on spatial patterns of rainfall. In further work, we plan to develop a more general description that will incorporate representation of the temporal and spatial rainfall structures (e.g., [55,57]). Such descriptions should include: a rain cell tracking algorithm (such as [23]), conceptual modeling of rain cell movement, and also changes in characteristics. In addition, the cell model can be extended to represent more complex structures such as ellipsoid or a more general decay function (e.g., [3,40,13,59,53]).