The hydroclimatology, hydrometeorology, and hydrology of flash floods in the arid/semiarid southwestern United States are examined through empirical analyses of long-term, high-resolution rainfall and stream gauging observations, together with hydrological modeling analyses of the 19 August 2014 storm based on the Kinematic Runoff and Erosion Model (KINEROS2). The analyses presented here are centered on identifying the structure and evolution of flood-producing storms, as well as the interactions of space–time rainfall variability and basin characteristics in determining the upper-tail properties of rainfall and flood magnitudes over this region. This study focuses on four watersheds in Maricopa County, Arizona, with contrasting geomorphological properties. Flash floods over central Arizona are concentrated in both time and space, reflecting controls of the North American monsoon and complex terrain. Thunderstorm systems during the North American monsoon, as represented by the 19 August 2014 storm, are the dominant flood agents that determine the upper tail of flood frequency over central Arizona and that also shape the envelope curve of floods for watersheds smaller than 250 km2. Flood response for the 19 August 2014 storm is associated with storm elements of comparable spatial extent to the drainage area and slow movement for the three compact, headwater watersheds. Flood response for the elongated and relatively flat Skunk Creek highlights the importance of the spatial distribution of rainfall for transmission losses in arid/semiarid watersheds.
Arid/semiarid regions, covering a large fraction of the global land area (Zeng et al. 2008), are affected by high-impact flooding associated with infrequent extreme rainfall events (e.g., Schick 1988). The hydrology, hydrometeorology, and hydroclimatology of flooding in arid/semiarid regions are poorly understood (e.g., Higgins et al. 2003), and yet are of great societal importance.
Our study focuses on the central Arizona region in the southwestern United States. Much of this region experiences 40%–60% of the mean annual rainfall during summer under the influence of the North American monsoon (NAM; e.g., Douglas et al. 1993). Flooding in the southwestern United States has been linked to warm-season convective rainfall associated with the North American monsoon (e.g., Gochis et al. 2006; Saharia et al. 2017), large-scale winter storms (Ely et al. 1994; House and Hirschboeck 1997; Neiman et al. 2013), and tropical cyclones (principally in the southwestern portion of the region). The relative importance of these different flood-generating agents in the southwestern United States varies with basin scale and geographical location (Gochis et al. 2006). In this study, we will focus on extreme flooding from convective rainfall, principally from warm-season organized thunderstorm systems during the NAM. These events are important flood agents for many watersheds in our study region, and they have posed particular challenges in examining flood response, because of the large variability of rainfall in time and space (Syed et al. 2003; Goodrich et al. 1995; Morin et al. 2006, 2007).
On 19 August 2014, multiple waves of thunderstorms crossed south-central Arizona and dropped heavy rainfall during a 12-h period leading to widespread flash flooding. The 19 August 2014 storm was characterized by strong forcing from an approaching upper-level trough over an unstable and moist atmosphere across central Arizona (FCDMC 2014). The upper-level trough extended toward the south off the coast of central California, with the contour of 5760 gpm at 500 hPa cut off from the main circulation and slowly evolving toward the east (FCDMC 2014). “Cutoff low” systems have been responsible for extreme rainfall and flooding in this region and other regions around the world (e.g., Knippertz and Martin 2007; Marquardt Collow et al. 2016; Awan and Formayer 2017; Pinheiro et al. 2016). The strong synoptic forcing served as a catalyst to initiate thunderstorms during the early morning hours on 19 August 2014. The atmosphere over Phoenix was humid and unstable prior to the initiation of convection. Precipitable water ranged from 47 to 52 mm through the early morning hours, which is approximately 1.75 times higher than the climatological mean for the period. The atmospheric profile was characterized by high convective available potential energy of 1400 J kg−1 and low convection inhibition of 33 J kg−1 at 1200 UTC (0500 LST) 19 August 2014 (FCDMC 2014). This storm was one of the most extreme rainfall events in recent decades, and it dominates the upper tail of the extreme rainfall and flood frequency distribution over Arizona. The 19 August 2014 event will be examined as a prototype for organized thunderstorm systems that are responsible for extreme flooding over the arid/semiarid regions of the southwestern United States.
Complexities of flood response over arid/semiarid watersheds are closely tied to extreme variability in rainfall and heterogeneities within the watersheds (e.g., vegetation, drainage networks, morphology, and bedrock type and state). Previous studies also identified the role of interactions between space–time rainfall variability and basin characteristics in determining scale-dependent flood response, for instance, storm size and motion relative to the drainage area and orientation of the drainage network (e.g., Smith et al. 2002; Borga et al. 2007; Yakir and Morin 2011; Morin et al. 2006, 2007; Morin and Yakir 2014). For arid/semiarid regions, storm size is typically small, and thus only part of the large watersheds can contribute to even the most extreme events. Transmission loss (i.e., reduction of streamflow through surface–aquifer interactions; e.g., Walters 1990; Vivoni et al. 2006) complicates the characterization of flood response in arid/semiarid watersheds by further increasing its dependency on temporal and spatial rainfall variability. The relative roles of rainfall variability and basin characteristics in determining flood response in arid/semiarid watersheds need further examination.
We address the following research questions in this study:
What are the characteristic patterns of storm structure and evolution for extreme flood-producing storms in arid/semiarid regions?
How does extreme flood response in arid/semiarid watersheds depend on temporal and spatial rainfall variability and watershed properties?
What is the nature of the upper tail of rainfall and flood magnitudes in arid/semiarid watersheds?
We examine these questions based on empirical analyses of observations from the Weather Surveillance Radar-1988 Doppler (WSR-88D) radar at Phoenix, high-resolution (15 min) rain gauge observations, and direct/indirect stream gauging records. Storm-tracking algorithms will be used to characterize the Lagrangian properties of flood-producing storms. We further examine the flood response of the 19 August 2014 storm based on hydrological modeling analyses using the Kinematic Runoff and Erosion Model (KINEROS2; Woolhiser et al. 1990; Goodrich et al. 2012).
We focus on four watersheds (each less than 200 km2, see Fig. 1b for locations) over Maricopa County, Arizona. The four watersheds exhibit contrasting geomorphological properties (e.g., mean elevation, slope, and drainage density; see Table 1 for details). Skunk Creek, New River, and Cave Creek are comparable in basin scales. The Cline Creek watershed is a subwatershed of Skunk Creek and is substantially smaller than the other three watersheds. Cline Creek, Cave Creek, and New River are headwater watersheds, with relatively large topographic gradients and compact shape. Skunk Creek, on the other hand, is flatter and elongated, with a larger drainage density than its counterparts (i.e., Cave Creek and New River). Skunk Creek has a higher density of residential areas; otherwise, the four watersheds exhibit little heterogeneity in land use patterns. There are contrasting spatial patterns of soil properties (i.e., texture and hydraulic properties) across the four watersheds (Fig. 1d, Table 2). The four watersheds are distributed over a region with large topographic gradients that experiences relatively frequent extreme rainfall (see analyses in section 3). Floods generated over the four watersheds pose significant social and economic impacts for the Phoenix metropolitan region that is located immediately downstream.
The paper is structured as follows. In section 2, we describe streamflow and rainfall observations, along with methods for radar rainfall estimation and storm tracking. Rainfall metrics representing flood response and the KINEROS2 model are also introduced in this section. The hydroclimatology of flash flooding is presented in section 3, followed by analyses of hydrometeorology (section 4) and hydrologic response (section 5) for the 19 August 2014 storm. A summary and conclusions are presented in section 6.
2. Data and methods
a. Streamflow and rainfall observations
The Flood Control District of Maricopa County (FCDMC) maintains an exceptional network of hydrometeorological observations, making this region an ideal setting for examining flood hydrology in arid/semiarid regions. There are 311 rain gauges distributed over the study region (see Fig. 1a for locations), with continuous records of 15-min rainfall accumulations from the early 1990s to the present. FCDMC provides stage data and discharge based on locally developed “stage discharge” rating curves for several major stream channels and tributaries.
Postflood surveys are carried out when the gauging stations are damaged or destroyed during extreme flood events. For instance, the stream gauge at the outlet of New River (Fig. 1b) was swept away near the time of peak runoff, and the channel was scoured by almost 3 m during the 19 August 2014 flood (K. Fossum, USGS, 2017, personal communication). A postflood indirect discharge measurement estimated a peak discharge of 923 m3 s−1 for this watershed. The U.S. Geological Survey also conducted an indirect discharge measurement over Cave Creek during the 19 August 2014 flood and reported a peak discharge of 286 m3 s−1. Continuous discharge observations and postflood records form the basis for the hydrological analyses of flood response.
FCDMC archives detailed statistics for each runoff event across the region, including the starting and ending times, duration, discharge peak timing, and magnitude. We term these “runoff events” because there is often an extended period of very low flow or no flow, common in intermittent and ephemeral tributaries in the region, between large precipitation events that trigger the recording of the gauging instruments. There are approximately 7900 runoff events across all the stream gauging stations. The record length is uneven across the stream gauging stations: 80% of the stations have records of more than 10 years, while 43% of the stations have records of more than 20 years. These long-record stream gauging stations provide a good representation of the topography and basin scale across the study region.
b. Radar rainfall estimation
Radar rainfall fields are derived using volume scan reflectivity observations from the WSR-88D radar located in Phoenix, Arizona (KIWA; see Fig. 1 for location). The KIWA radar suffers from beam blockage by terrain in the northeast sector that affects the first two tilts (0.5° and 0.9°). There is only a narrow band of blockage to the north of the radar for the second tilt (0.9°). We confine our study region over Maricopa County and use the second tilt for rainfall estimation. The third tilt (1.4°) is used where beam blockage occurs for the second tilt. We do not use the first tilt because of beam blockage and also contamination by ground clutter over the study region (figure not shown).
Conversion of reflectivity to rainfall rate is based on Z–R relationships taking the form , where R is rain rate (mm h−1) and Z is the radar reflectivity factor (mm6 m−3). The choice of the Z–R parameters can impose significant bias on radar rainfall estimates. We use the standard NWS Z–R relationship (a = 0.017, b = 0.714), as well as a cap of 53 dBZ to reduce unreasonably large estimates caused by hail and graupel cores in thunderstorms (Fulton et al. 1998). We use NOAA’s Weather and Climate Toolkit (https://www.ncdc.noaa.gov/wct/) to convert rainfall fields from a 2D polar coordinate system to a Cartesian coordinate system with horizontal resolution of approximately 1 km. The converted rainfall fields are in irregular temporal resolutions of 5–6 min. We linearly interpolate rainfall fields to 1 min and then integrate into hourly and daily scale for future applications.
We implement daily mean field bias correction (Smith and Krajewski 1991; Krajewski and Smith 2002) using the network of 311 rain gauges to improve the rainfall estimates by reflectivity-based radar. The daily bias correction takes the form
where is the daily rainfall accumulation for gauge j on day i, is the daily rainfall accumulation for the radar pixel containing gauge j on day i, and is the index of the rain gauge stations where both rain gauge and radar have positive rainfall accumulation for day i. The bias value is set to be different from 1.0 only if there are at least 10 positive radar–rain gauge pairs in the domain. Bias correction is implemented by multiplying each 15-min rainfall field on day i by . The daily bias correction removes systematic bias due to variability in Z–R relationships and radar calibration errors (Villarini and Krajewski 2010). This bias correction procedure is the same as that used in, for example, Wright et al. (2014), Yang et al. (2013), and Smith et al. (2012).
The bias correction scheme was implemented over 118 warm-season days (July–September) spanning 1995–2014. The storm days were determined based on rain gauge observations and the completeness of radar coverage. A “day” is defined from 1200 to 1200 UTC of the next day. The daily correlation coefficient (CC) improved from 0.7 to 0.81 with bias correction (Fig. 2). The daily root-mean-square error (RMSE) decreased from 10.0 to 7.0 mm, and the daily mean absolute error (MB) decreased from 5.0 to 3.0 mm. A scatterplot between the rainfall accumulations from the 311 rain gauges and corresponding radar-pixel rainfall estimates for the 19 August 2014 storm illustrates the error characteristics of rainfall estimates (Fig. 3). The correlation coefficient for the hourly accumulation is 0.67; for the 12-h accumulation, the correlation coefficient is 0.83. The maximum 12-h rainfall accumulation of 137 mm was captured by the bias-corrected radar rainfall fields.
The improvements in radar rainfall estimates with multiplicative bias correction are consistent with previous studies showing that mean field bias correction is a useful step in developing accurate rainfall products from conventional radar reflectivity measurements (Krajewski and Smith 2002; Gourley and Vieux 2005; Smith et al. 2012). We note that comparisons between rain gauge and corresponding radar-pixel rainfall estimates at a 1-h time scale still contain significant errors (Fig. 3a). Detection of the error sources can be difficult, even though we believe errors could be partially due to the use of the conventional Z–R conversion schemes. Upgrade of the U.S. Next Generation Weather Radar (NEXRAD) network to dual-polarization capabilities provides an alternative for radar rainfall estimates using a more sophisticated algorithm (see, e.g., Ryzhkov et al. 2005; Gourley et al. 2010). Performance of radar rainfall estimates in Phoenix based on dual-polarization products will be evaluated in future studies.
c. Storm-tracking algorithms
We utilize the Thunderstorm Identification, Tracking, Analysis, and Nowcasting (TITAN; Dixon and Wiener 1993) storm-tracking algorithms to examine key aspects of storm structure and evolution for the 19 August 2014 event. Storm-tracking analyses are based on 3D reflectivity fields obtained from the KIWA radar in Phoenix (see Fig. 1a for location). Time-varying properties for each storm element are derived as a part of the identification and tracking procedures. The variables analyzed in this study include the area of storm elements (projected area of the storm envelope when viewed from above), maximum reflectivity, storm speed, and storm direction (see, e.g., Yang et al. 2016; Yeung et al. 2015).
We use a 45-dBZ reflectivity threshold and 5-km3 volume threshold to identify storm elements for each storm event. A contiguous region of reflectivity values exceeding 45 dBZ is identified as a storm cell provided that the total volume with reflectivity exceeding 45 dBZ is also greater than 5 km3. Echo-top height for each storm element is defined as the maximum height of radar reflectivity exceeding 18 dBZ. We use the term “storm element” for the features tracked using the above thresholds. Each storm element is automatically fitted by an ellipse providing a set of geometrical parameters, including geographic location of the centroid, major/minor axis, and orientation of the ellipse. Storm elements responsible for rainfall and flood over the four watersheds are examined in detail.
d. Rainfall-weighted flow distance
We analyze flood response of the 19 August 2014 storm based on time series of rainfall properties that are linked to drainage network structure (see, e.g., Smith et al. 2002). The drainage network for each watershed is represented by the distance function . The value of for each point x within the drainage basin is computed as the sum of the overland flow distance from x to the nearest channel and the distance along the channel to the basin outlet.
The normalized flow distance time series is a function of the rainfall field and the flow distance function (Smith et al. 2002). It is defined as the ratio of the rainfall-weighted centroid distance to the basin outlet and the maximum distance from the basin outlet (including flow path from both overland and channel routing). The distance time series takes the form
where A is the spatial domain of the watershed, and the weight function is given by
Values of range from 0 to 1, with values close to 0 indicating that rainfall is distributed near the outlet and values close to 1 indicating the far periphery of the basin. If rainfall is uniformly distributed over the basin, then the weights do not vary spatially, and we obtain the mean flow distance:
The rainfall-weighted flow distance dispersion (Smith et al. 2002) is defined by
where is the dispersion for uniform rainfall:
The normalized dispersion takes the value 1 for uniform rainfall, is less than 1 for rainfall characterized by a unimodal peak, and can take values greater than 1 (e.g., in the case of multimodal rainfall peaks in both the upper and lower basin).
We further examine flood response of the 19 August 2014 storm over the four watersheds based on KINEROS2. KINEROS2 is an event-oriented, physically based rainfall–runoff model developed for watersheds in arid/semiarid environments (http://www.tucson.ars.ag.gov/kineros/; Woolhiser et al. 1990; Goodrich et al. 2012). The watershed is represented by a cascade of overland flow model elements (planar or curvilinear) and channel model elements, in which flow and sediment are routed from one overland element to another and ultimately to the channel. KINEROS2 uses the Parlange three-parameter model to estimate infiltration (Smith and Parlange 1978). Infiltration may occur from either rainfall directly on the soil or from the ponded surface water created from upslope rainfall excess. Infiltration can also occur over the channel bed during the routing processes, representing the reduction of streamflow through surface–aquifer interactions (i.e., transmission loss). Routing over the planes or in the channels is described by the one-dimensional kinematic wave equation. The numerical solution provides discharge at any point in time and at any distance along a flow path.
The Agricultural Research Service of the U.S. Department of Agriculture (USDA-ARS), the University of Arizona, and the U.S. Environmental Protection Agency developed an Automated Geospatial Watershed Assessment (AGWA) tool that can be used to generate parameters for KINEROS2 based on GIS data layers of topography, soil, and land use/land cover (https://www.tucson.ars.ag.gov/agwa/; Goodrich et al. 1997; Miller et al. 2007). Setup of KINEROS2 in this study was based on the 1-arc-s (approximately 30 m) National Elevation Dataset (available at http://nationalmap.gov/) for watershed delineation and stream location. We used the State Soil Geographic database (STATSGO, available at https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/survey/geo/?cid=nrcs142p2_053629) for soil data and the 2011 National Land Cover Database (https://www.mrlc.gov/nlcd2011.php) for land use/land cover in the model (see Fig. 1 and Table 2 for details). Model calibration and validation will be introduced in section 5b.
3. Hydroclimatology of flash flooding
The 7900 runoff events exhibit a “bimodal” seasonal distribution, with a relatively larger peak (a total number of 3698) occurring in summer and early autumn (i.e., July–September, referred to as NAM events) and a second peak (a total number of 3038) occurring during the period from winter to early spring (i.e., December–March, referred to as winter events; Fig. 4). Approximately 20% of the total events are observed in August.
The distributions of peak magnitudes between NAM and winter flood events are different. There is a heavier tail toward extreme events of large flood peak magnitudes (exceeding 0.1 m3 s−1 km−2) during NAM, while more frequent events with smaller magnitudes (less than 0.01 m3 s−1 km−2) occurred in winter (Fig. 5). Figures 4 and 5 reflect the dominant role of thunderstorm systems during the North American monsoon as the major flood agents over central Arizona.
We examine the dependence of peak unit discharge (i.e., peak discharge divided by drainage area) on drainage area. As can be seen from Fig. 5, there is a pronounced decreasing trend of peak unit discharge with drainage area, especially for large peak magnitudes, that is, the “envelope curve” portion of the plot. The decreasing trend is mainly due to spatial and temporal averaging of the flood-generating rainfall intensities with increasing drainage areas. The reduced peak unit discharge is also partially due to the losses of flood water into channel alluvium, especially for large watersheds with relatively long flow distances. As noted by Goodrich et al. (1997), in these influent hydrologic environments, runoff response per unit area decreases with increasing drainage area and becomes more nonlinear based on analysis of the Walnut Gulch Experimental Watershed data in southeastern Arizona. Similar features were also identified in the eastern Mediterranean region (e.g., Borga and Morin 2014).
Another feature pertaining to the envelope curve is the shift of dominant controls on the shape with the increase of drainage area (Fig. 5). For small watersheds (smaller than 250 km2), the envelope curve is predominantly controlled by NAM flood events. There are 129 NAM flood events with peak magnitudes exceeding 1 m3 s−1 km−2 over small watersheds, while only six winter flood events fall into this scale range. The NAM events are primarily caused by thunderstorms of limited spatial extent. This factor combined with channel transmission losses favors higher peak unit discharges for smaller watersheds (Goodrich et al. 1997; Syed et al. 2003). The shape of the envelope curve of peak unit discharge for large watersheds (larger than 250 km2) reflects a mixture of both winter and NAM flood events, with relatively more peaks occurring in winter. Winter floods are frequently linked to frontal systems with high rainfall totals but lower intensities occurring over larger areas, melting of snow, and/or rain on snow (House and Hirschboeck 1997). The contrasting controls of NAM and winter flood events highlight the dependence of flood response on the intrinsic properties of flood-producing storms and drainage basins and will be analyzed in the following sections.
Flood events during the North American monsoon in central Arizona are concentrated in both space and time. The spatial distribution of stream gauging stations (Fig. 6) that recorded the 129 NAM flood events with peak magnitudes exceeding 1 m3 s−1 km−2 (see also Smith and Smith 2015) highlights the spatial distribution of floods. The events are distributed over 37 watersheds, with drainage areas ranging from 2 to 250 km2 (average size is 49 km2). Most of the 37 watersheds are located within or close to the Phoenix metropolitan region and extend toward the elevated terrain north of the city. The 129 flood events are produced by 53 individual storm episodes, resulting in approximately 2–3 flood events per storm episode. An individual storm episode is defined as a period of consecutive rainy days separated from other storm episodes by at least two nonrainy days. A “rainy day” is determined based on rainfall observations from all 311 rain gauges together with radar observations. There are specific storm systems that are responsible for a large number of extreme flood events. The 19 August 2014 storm, for example, produced 18 flood events with peak magnitudes exceeding 1 m3 s−1 km−2 (Fig. 5).
The mean daily rainfall for 25 storm days over central Arizona illustrates the spatial heterogeneity of flood-producing rainfall (Fig. 7). The 25 days are selected from the catalog of 53 storm episodes (due to the availability of radar rainfall fields). Some of the most severe storm episodes (e.g., the 19 August 2014 storm) are included in the 25 storm days. The maximum rainfall (as represented by the 20-mm contour in Fig. 7) is located over the mountainous region that includes the headwaters of New River and Cave Creek, and the 18-mm contour extends toward the southeast, following the sharp gradients of the topographic transition (see also Fig. 1b for high-resolution local topography). The axis of maximum daily rainfall reflects the role of topography in determining the spatial variability of flood-producing storms. Orographic thunderstorm systems are important elements of extreme flood climatology (e.g., Laing and Fritsch 1997; Javier et al. 2007). A recent study by Mascaro (2017) shows that rainfall maxima over central Arizona are controlled by elevation for all durations during the monsoon season, and summer thunderstorms are localized in both space and time. The spatial pattern of mean daily rainfall for the 25 storm days obtained in this study is consistent with the spatial distribution of climatological mean daily rainfall over central Arizona (see Fig. 7b in Mascaro 2017).
4. Hydrometeorology of the 19 August 2014 storm
a. Rainfall statistics
Heavy rainfall for the 19 August 2014 event extended across a 12-h period from 1200 UTC August 19 to 0000 UTC 20 August 2014. The storm total rainfall accumulations based on both rain gauges and radar rainfall estimates (Fig. 8) illustrate the concentration of rainfall over complex terrain. The maximum gauge rainfall accumulation is 131 mm (gauge labeled as “A” in Fig. 8), while the maximum radar rainfall accumulation is 165 mm. Three of the four gauges with total accumulation exceeding 100 mm are located in the New River basin. The fourth is near the outlet of Cline Creek (labeled as “D” in Fig. 8). Detailed rainfall statistics for the four rain gauges with the most extreme rain rates are summarized in Table 3. The estimated recurrence intervals for the maximum 6-h rainfall accumulation exceed 100 years for all four gauges based on the NOAA Atlas 14 point frequency estimates. Gauge A experienced the most extreme rain rates at 6- and 12-h time intervals during 19 August 2014, with estimated recurrence intervals of 833 and 541 years, respectively. The rain gauge at the outlet of New River (labeled as “C” in Fig. 8) received 120 mm of rain in 6 h, with the estimated recurrence interval of 896 years.
b. Storm structure and evolution
There are approximately three storm periods within the 12-h heavy rainfall period on 19 August 2014. The first storm initiated around 1300 UTC and moved northeast across the four watersheds during the following 90 min. The second storm started right after 1500 UTC; small storm elements initiated along the western boundary of New River and moved toward the northeast. The third storm was characterized by a cluster of small storm cells and started around 2200 UTC. The first and second storm periods were responsible for extreme flooding over the four watersheds and will be analyzed below.
The time series of basin average rain rate derived from the radar-estimated rainfall fields illustrate rainfall forcing of floods in the four watersheds (Fig. 9a). The two storm periods exhibit contrasting features of storm evolution and structure. Storm elements in the first storm period started over New River and Skunk Creek and then moved toward the east (1326 UTC, Fig. 10a). There were also several small storm elements located in the southeast corner of Skunk Creek (1342 UTC, Fig. 10b). A storm merger occurred at 1355 UTC (Fig. 10c), leading to a larger storm cluster with a size of 160 km2, which is comparable to the scales of the three large watersheds (Fig. 9c). The large storm element moved across Cline Creek, New River, and Cave Creek successively during the following hour. Intensification of the storm core over New River can also be noted during its motion toward the mountainous terrain (1411 UTC, Fig. 10d), even though the maximum radar reflectivity for the storm elements during the first storm period does not show significant increasing trends and maintains an almost constant value around 55 dBZ (Fig. 9b).
For the second storm period, storm elements initiated to the west of New River and moved toward the northeast (80°, zero is toward true north; see Figs. 9d and 11a). Storm cells repeatedly initiated and persistently passed through New River, exhibiting features of training storm systems (characterized with a repeated motion of cells over a particular region; see, e.g., Doswell et al. 1996). The training storm elements were associated with a strong steering-level wind (exceeding 13 m s−1) and terrain blocking. Storm cell size and intensity were enhanced during the propagation toward the mountainous terrain (Figs. 9b,c). The strong convergence zone over the windward slope is favorable for the initiation of convective cells. Preferential initiation of storms in specific locations and upscale growth of convective systems are common features of storm structure and motion in complex terrain (Smith et al. 1996; Laing and Fritsch 1997; Javier et al. 2007; Tarolli et al. 2013). Another notable feature of the second storm period is the association of peak surface rain rates with the decreased storm speed. The storm speed decreased from 40 to less than 25 km h−1 within 34 min (1624–1658 UTC; see Figs. 9d and 11a–c). The slower motion of storm elements produced a steady peak rain rate over New River and was responsible for 6-h rainfall accumulation exceeding 120 mm.
Key properties of storm elements that passed through each of the four watersheds for the 19 August 2014 storm are characterized by storm scale, convective intensity, and storm motion (Fig. 12). The median values of storm cell sizes are around 10 km2 for the four watersheds, indicating that there are numerous “pop-up” storm cells during the event. The mean values of storm cell size (around 40 km2) are comparable across the three large watersheds (Skunk Creek, New River, and Cave Creek). Mean storm cell size for Cline Creek is 60 km2, more than twice as large as the basin scale (28 km2). Maximum reflectivity ranges from 50 to 55 dBZ (with the mean values around 52 dBZ, corresponding to 88 mm h−1 using the current Z–R relationships) and varies marginally across the four watersheds. The median speed of the storm cells is around 25 km h−1, with a dominant direction of motion toward the northeast (around 70°). The storm cell size is comparable with that derived from a climatological study over northern Israel (using a different threshold of 40 dBZ), with mean values ranging from 50 to 64 km2 (Peleg and Morin 2012), though the maximum reflectivity exhibits comparable ranges.
5. Flood response of the 19 August 2014 storm
a. Rainfall variability and drainage network
Figure 13 shows the observed and simulated hydrographs (based on KINEROS2 and which will be elaborated on in the following subsection) of the 19 August 2014 flood event for the four watersheds. There are two separate flood peaks over the four watersheds associated with the distinct storm episodes (see the yellow dots in Fig. 13). For New River and Cave Creek, the second flood peaks are larger than the first, while the opposite is true for Skunk Creek and Cline Creek. The relative magnitudes of flood peaks are related to the storm motion, volume, and/or spatial extent of storm core (e.g., Syed et al. 2003; Morin and Yakir 2014). The consistency between the time series of basin-averaged rain rate and spatial coverage of heavy rainfall is notable (Figs. 14a,b), indicating strong correlations between the spatial extent of the storm core and basin-averaged rain rate. The heavy rainfall area covered almost the entire Cline Creek region and produced a peak unit discharge of 5.8 m3 s−1 km−2. The spatial coverage of heavy rainfall over Skunk Creek, less than 30%, is the smallest (Fig. 14a), leading to smaller maximum 15-min, maximum 60-min, and maximum 2-h rain rates compared to the other three watersheds (Table 4).
We focus on storm motion over a 2-h time window that experienced a maximum 2-h basin-averaged rain rate during the two storm periods (see the shaded area in Fig. 14). The time derivative of implies storm motion relative to the drainage network, with positive (negative) values indicating upstream (downstream) motion. The first flood peaks are associated with storm elements moving upstream relative to the drainage networks, especially for the three headwater watersheds [New River, Cave Creek, and Cline Creek, as increases from 1300 to 1500 UTC; see Fig. 14]. The second flood peaks are associated with the slow motion of storm elements relative to the drainage networks [or slightly downstream motion for Cave Creek, as does not show notable changes with time or a slightly decreasing trend for Cave Creek; Fig. 14]. Slow-moving storm elements are key features for extreme flood response (e.g., Doswell et al. 1996; Morin and Yakir 2014; Yang et al. 2016).
Peak basin-averaged rain rates are associated with uniformly distributed rainfall fields for the three headwater watersheds, while the storm core for Skunk Creek is relatively smaller than the basin scale [with the dispersion value smaller than 1.0; see Fig. 14a] and predominantly centered over the upstream region during the event [with the normalized flow distance larger than the mean value ]. The normalized flow distance for storm total rainfall is close to the mean flow distance (with a uniform rainfall field) over the three headwater watersheds (Table 4), indicating the strong smoothing effects of drainage networks on spatial rainfall variability in producing extreme floods. For Skunk Creek, the smoothing effect is less significant because of the large deviation between the normalized flow distance for storm total rainfall (0.72) and mean flow distance with uniform rainfall distribution (0.59; Table 4). The first flood peak over Skunk Creek, with peak unit discharge of 2.3 m3 s−1 km−2, is related to the propagation of the flood hydrograph from upstream to the outlet, exhibiting a steep rising limb and larger time delay between the flood peak and rainfall peaks (approximately 3 h, as opposed to less than 1 h for the other three watersheds; Fig. 13).
b. Model calibration and validation
We calibrated KINEROS2 focusing on Manning’s roughness coefficient and saturated hydraulic conductivity for both the channel and hillslope elements. The other parameters (e.g., geometric parameters and soil hydraulic properties) were derived from the GIS data layers (topography, soil, and land use) using AGWA. The saturated hydraulic conductivity for the channel elements was set to 80 mm h−1 (Rawls et al. 1982; Rawls and Brakensiek 1985). The saturated hydraulic conductivity for hillslope elements was automatically determined using AGWA by averaging values over different soil types within the element. We further adjusted both Manning’s roughness coefficient and the saturated hydraulic conductivity by applying multipliers for each of the two parameters to all the elements in the model while maintaining the original heterogeneity across the watersheds. The model was calibrated over Skunk Creek for the 19 August 2014 flood event, with the calibrated multipliers directly transferred to the other three watersheds for model validation. The initial wetness for the watershed was set to 0.098, representing relatively dry conditions over the entire region before the storm.
The calibration is based on the Parameter Estimation Tool (PEST; Doherty et al. 1994), with the objective function to minimize the sum of squared errors between the simulated and observed hydrographs at the basin outlet. The calibrated multiplier for the saturated hydraulic conductivity and Manning roughness coefficient are 0.78 and 1.4, respectively. Figure 13a evaluates the performance of the model over Skunk Creek for the 19 August 2014 flood. The Nash–Sutcliffe efficiency coefficient (Nash and Sutcliffe 1970) is 0.68, indicating a good consistency between the simulated and observed hydrograph of the 19 August 2014 flood.
We further evaluate the performance of the model with calibrated multipliers over the other three watersheds. Evaluation statistics are summarized in Table 5. KINEROS2 accurately captured the flood peak timings of this event over each of the four watersheds, with errors ranging from 0 (Cline Creek) to 0.6 h (Skunk Creek). The large error in Skunk Creek is partially related to the long response time (about 3 h, i.e., time delay between the rainfall peak and flood peak). There are additional discrepancies between the simulated and estimated peak discharges. For instance, the relative biases for Cave Creek and Skunk Creek are 34% and 24%, respectively. However, when the uncertainty range of the indirect peak discharge measurements is considered, the simulated flood peak over Skunk Creek falls within the uncertainty range, and the simulated flood peak over Cave Creek falls within 27 m3 s−1 of the upper range of uncertainty. We note, however, that the model reproduced flood peak magnitudes reasonably well over New River and Cline Creek, with relative biases of less than 15%. If the quality of the FCDMC indirect measurements were comparable to that of the USGS indirect measurements over Skunk Creek and Cave Creek, the simulated peaks for New River and Cline Creek would fall well within the uncertainty range.
Simulation of extreme flood processes over arid/semiarid regions is challenging (e.g., Michaud and Sorooshian 1994; Pilgrim et al. 1988; Vivoni et al. 2006) because of a poor understanding of the physical mechanisms of runoff generation and flood routing, uncertainties resulting from errors in radar rainfall estimates, uncertainties in model parameters, etc. Because of the losing or influent character of arid/semiarid watersheds, the runoff-to-rainfall (signal to input) ratio can be very small. Thus, uncertainties in the rainfall inputs (noise) can be a very large percentage of the runoff signal (Yatheendradas et al. 2008; Goodrich et al. 2012). Another challenge stems from observations of flash floods in arid/semiarid regions (e.g., Borga et al. 2010, 2014; Berg et al. 2016; Halbert et al. 2016): debris flow, erosion of the river channel, scour on the rising limb of the hydrograph and deposition on the falling limb (rating curves are thus subject to change from event to event and dynamically within an event), and damages to the gauging stations can all adversely affect streamflow measurements. As noted, the New River stream gauge was destroyed by the 19 August 2014 flood near the time of peak runoff. The remaining portion of the hydrograph was estimated using recession curves common to the streams in the area (K. Fossum, USGS, 2017, personal communication). Therefore, less weight should be given to the model underprediction of the recession limb of this storm hydrograph (Fig. 13b).
Despite these uncertainties, KINEROS2 captured key features of flood response (in terms of peak magnitude and peak timing) for the 19 August 2014 storm, with good consistency between the simulated and observed hydrographs over the four watersheds (correlation coefficients around 0.80; Table 5). We will rely on the model to analyze the hydrological response of the flood event.
c. Effects of rainfall variability and distributed watershed properties
Figure 15 summarizes the storm-event water balance for the 19 August 2014 flood over the four watersheds based on the model simulations. The storm-event runoff ratio (i.e., total outflow volume to total rainfall accumulation) is the largest for New River (0.42), with Cave Creek being the smallest (0.22). The small runoff ratio for Cave Creek is mainly attributed to the largest infiltration ratio (i.e., total infiltration volume to total rainfall accumulation) over the plane elements (0.70). A notable feature is that Skunk Creek has the largest infiltration ratio (0.20) over the channels, more than twice as large as that of the other three watersheds (0.06–0.08), indicating a significant contribution of transmission loss in determining the storm-event water balance over Skunk Creek for the 19 August 2014 flood.
Figure 16 shows the spatial distribution of runoff ratios and infiltration ratios of the 19 August 2014 flood for each plane and channel element over the four watersheds, exhibiting the roles of rainfall variability and distributed watershed properties (i.e., soil textures, land use/land cover, and bedrock) in determining the water balance and then flood peak response. Plane elements with storm-event runoff ratios exceeding 0.80 are distributed over the upstream portion of Skunk Creek, the western portion of New River, and the middle to lower portion of Cline Creek (Figs. 16a,c,d), which is consistent with the spatial pattern of storm total rainfall accumulation over each of the three watersheds (Fig. 8). Larger runoff ratios can be partially related to the role of distributed soil hydraulic properties over the watersheds in dominating runoff generation processes. For instance, the upstream portion of Skunk Creek is covered with soils composed of gravelly units and unweathered bedrock that produce a low saturated hydraulic conductivity, large rock fraction, and porosity (also similar to New River; Fig. 1d, Table 2). Cave Creek, which is mainly covered by loam of high saturated hydraulic conductivity, exhibits relatively smaller runoff ratios (<0.60). For Skunk Creek, the main runoff production is concentrated over the upstream region, leading to a longer routing time for the runoff to reach the outlet than the other three watersheds. Long response time facilitates the loss of storm flow into river channels through surface–aquifer interactions (i.e., transmission loss). The reinfiltrated runoff through channel beds takes up more than 15% of the total runoff generation over Skunk Creek, especially through the channels at the downstream region where runoff ratios are smaller than 0.20 (Fig. 16c), indicating the important role of transmission loss in dictating flood response over Skunk Creek. There are some channel elements exhibiting infiltration ratios close to 100% over Skunk Creek. The infiltration ratios over the channels are relatively small over New River, Cave Creek, and Cline Creek, that is, making up less than 10% of total runoff over most of the channel elements. Plane elements with large runoff ratios are distributed close to the outlets of New River and Cave Creek, paired with the spatial distribution of channels of low infiltration ratios. Contrasting spatial distributions of runoff ratios and infiltration ratios for the storm-event water balance highlight the importance of improved representations of spatial rainfall variability and physical processes for transmission loss in characterizing flood response in arid/semiarid watersheds.
6. Summary and conclusions
The hydroclimatology, hydrometeorology, and hydrology of flash flooding in the arid/semiarid southwestern United States are examined through empirical analyses of rainfall and stream gauging observations together with hydrological modeling analyses of an extreme storm event that produced widespread flash floods over the region. We focus on the properties of flood-producing storms that characterize the upper tail of rainfall and flood magnitudes in four arid/semiarid watersheds in Maricopa County, Arizona. The principal conclusions of this study are summarized as follows:
Flash floods over central Arizona are concentrated in both time and space. There is a distinct bimodal seasonal distribution of the runoff events over this region, with more than 40% of the events occurring during the North American monsoon. Flood events with peak magnitudes exceeding 1 m3 s−1 km−2 are concentrated over small watersheds (less than 250 km2) that are predominantly located within the Phoenix metropolitan region or the surrounding region toward the mountainous terrain to the north of the city. Several severe storm events during the North American monsoon season are responsible for a majority of extreme flash floods. Composite rainfall fields from those flood-producing storms exhibit pronounced spatial heterogeneity closely tied to strong orographic effects on precipitation.
Floods during the North American monsoon control the shape of the envelope curve for smaller drainage basins (less than 250 km2), while the shape of the envelope curve for larger drainage basins is dominated by floods during winter. Thunderstorm systems during the North American monsoon are the dominant flood agents that characterize the upper tail of extreme flood frequency distributions over central Arizona. Contrasting controls on the flood envelope curve reflect the role of interactions between the spatial scales of flood-producing storms and drainage basins in determining scale-dependent flood response.
The 19 August 2014 storm, which produced record rainfall accumulations and flood peak discharges, defines the upper-tail properties of extreme flood frequency over central Arizona. Extreme rainfall and flooding from this event was related to organized thunderstorm systems with strong forcing from an approaching upper-level trough imposed over an unstable and moist storm environment. Four rain gauges that recorded the most extreme 6- and 12-h rainfall accumulations are located over the topographic transition regions (i.e., windward regions), indicating that topography plays a critical role in initiating and enhancing convection through orographic uplifting.
Analyses based on storm-tracking algorithms are used to characterize key aspects of storm motion and structure for the 19 August 2014 storm. Convective storm elements frequently initiated to the west of the New River because of complex terrain and persistently propagated across the watershed, exhibiting typical features of training storm systems (e.g., Doswell et al. 1996). The slow motion of storm elements and the training effect were responsible for the extreme 6-h rainfall accumulations, with an estimated return period of 896 years over New River (gauge A in Fig. 8). Peak basin-averaged rain rates over the four watersheds were associated with storm elements with spatial coverage comparable to or larger than the drainage areas of the watersheds.
There are two flood peaks associated with two separate storm periods over each of the four watersheds. The first flood peaks were produced by storms moving upstream relative to the drainage networks, while the second flood peaks were associated with slow storm motion relative to the drainage networks. Skunk Creek exhibited a spatial concentration of heavy rainfall over the upstream portion of the watershed, leading to longer flow distances from the runoff-generation region to the outlet that enhanced transmission losses. A notable feature is the consistency of the fractional coverage of the heavy rainfall region over the four watersheds with basin-averaged rain rates, a key element in flood response for the three headwater watersheds.
Spatial distributions of runoff ratios and infiltration ratios for the storm-event water balance over the four watersheds highlight the important roles of rainfall variability and distributed watershed properties (i.e., soil textures, land use/land cover, bedrock) in determining the storm-event water balance and flood peak response. For Skunk Creek, the main storm flow is generated over the upper portion of the watershed due to the large rainfall accumulation and soil hydraulic properties of low infiltration intensity, and 20% of the generated storm flow reinfiltrates into the aquifers underlying the channel beds downstream during the propagation of the flood wave to the outlet. For the other three watersheds, the proximity of plane and channel elements with high runoff ratios and low infiltration ratios to the outlet limits the consumption of storm flow through transmission loss. Improved representations of spatial rainfall variability and physical processes for transmission loss are key elements in better characterizing flood response in arid/semiarid watersheds.
This research was jointly supported by the United States–Israel Binational Science Foundation (Grant BSF-2016953) and the National Science Foundation (Grants EAR-1632048, AGS-1522492, and CBET-1444758). The authors would like to acknowledge the Flood Control District of Maricopa County (FCDMC), Arizona, for maintaining the rainfall and stream gauging stations and making observational datasets publicly available.