The study area was limited to central Arizona to minimize geologic and climatic variation as well as differences in stormwater management regulations. This region includes the Phoenix metropolitan area. Flow records were obtained from USGS and Flood Control District of Maricopa County gages to capture a range of development, extent, and basin size. USGS data were obtained at 15-minute and daily intervals. Flood Control District data had variable time intervals, with higher-resolution data (e.g. 1 to 10 minute intervals) collected during periods of flow - these data were interpolated to 15-minute intervals and averaged at daily time steps. Sites that had at least 14 years of recent data were selected, and the data set was restricted to the 2003 to 2016 water years (where the water year was 1 October to 30 September). This number of years allowed for a sufficiently long period of time to capture climatic variability, while still having an adequate quantity of sites across a range of urban development. Gages that had more than a 10-cm gap documented between the channel bed and water-level sensor were excluded to most accurately characterize low flows. These criteria yielded 19 sites, which drained watersheds ranging in size from 1.5 to 425 km2.

Watersheds were characterized by landscape and climatic factors that could potentially explain observed relationships with hydrologic metrics. Per cent development and per cent impervious surface were based on the National Land Cover Database (Homer et al., 2015; Xian et al., 2011). Watershed development and impervious cover were calculated from National Land Cover Database rasters from 2001, 2006, 2011, and 2016, and for years in between. Total development was calculated by summing the high, medium, and low density development and developed open space categories. Watershed characteristics, such as area, annual precipitation, and slope were obtained from USGS StreamStats or Flood Control District gage information, mean area-weighted flowrate was calculated using streamflow records and watershed area. The total percentage of watershed area that was covered by stormwater control measures, such as retention basins, was calculated for a subset of watersheds with available data. Data regarding stormwater control measures are thought to be fairly accurate, though there are likely to be small errors in recorded basin area.

Various hydrologic metrics related to flow magnitude or flashiness were calculated for each watershed. The Indicators of Hydrologic Alteration software (The Nature Conservancy, 2009) was used to generate metrics based on daily flow data, since this is a software that has been commonly used to evaluate flow regime and ecological flow metrics (e.g., McDaniel and O’Donnell, 2019; Moltz et al., 2018; Schoonover et al., 2006; The Nature Conservancy, 2006). Both Indicators of Hydrologic Alteration and another increasingly used similar software—the USGS’s EflowStats package for R—only calculate metrics for daily flow data (US Geological Survey, 2017). However, it is also important to examine flow patterns at sub-daily intervals. In arid streams, flow events may last for less than a day for smaller storms and catchments; these would not be represented by daily data. Anthropogenic activities may also influence flows directly at certain times of day, such as via inputs of treated wastewater or irrigation-related runoff. Sub-daily and daily metrics for characterizing altered flow regimes can be different and not always even well-correlated (Bevelhimer et al., 2015), or may provide similar results (Schoonover et al., 2006). Therefore it is important to use sub-daily data particularly when trying to tease out impacts of human interventions (Bevelhimer et al., 2015). Sub-daily streamflow data available at selected sites allowed for a comparison using sub-daily and daily data.

Daily hydrologic metrics used in this study as output from Indicators of Hydrologic Alteration included zero flow days, rise rate, and fall rate. Zero flow days were calculated as an annual total of days with no measurable flow while rise and fall rates were calculated as an annual median rate. Rise and fall rate were calculated separately for sub-daily flow data using R (R Core Team, 2018). A similar algorithm as is used in the Indicators of Hydrologic Alteration software, where we examined the change in discharge between each timepoint and calculated the median annual rate of increase (rise rate) and decrease (fall rate), was applied. This algorithm was applied to the original (uninterpolated) sub-daily data to associate a given discharge change with the most accurate time interval. The Richards-Baker Flashiness Index (R-B Flashiness Index), coefficient of variation (CV), and peaks-over-threshold (POT) were calculated on both daily and sub-daily flow data using R, where R-B Flashiness and CV were calculated as annual values and POT was calculated as an annual sum of peaks over the selected threshold. R-B Flashiness Index and CV are both metrics for the variation in streamflow, where R-B Flashiness Index is calculated as the difference in volumetric flow rate between adjacent timepoints, normalized over the sum of volumetric flow rate for the entire year (Baker et al., 2004), and CV is the ratio of the standard deviation and mean volumetric flow rate for a given year. POT were calculated for flow data that were normalized by watershed area. There are multiple approaches to calculating thresholds for POT analysis. Specific numeric thresholds rather than percentiles derived from an aggregation of data from all sites were used given that a few watersheds are nested. Several thresholds were evaluated, with the maximum of 1 m3 s-1 km-2 chosen given that it was used in an analysis of stream flashiness across the U.S. using flow data at 15-minute intervals (Smith and Smith 2015). An additional lower threshold for sub-daily data (0.1 m3 s-1 km-2) was based on a target of 2 to 3 POTs per year per site on average (Mallakpour and Villarini, 2016; Mediero et al., 2014). Thresholds for POT analysis of daily data (0.01 and 0.05 m3 s-1 km-2) were chosen to have similar annual frequencies to the sub-daily thresholds. For sub-daily data, only one POT was counted per six-hour interval. For daily data, each day exceeding the chosen threshold was counted as a POT. Given the ephemeral nature of streamflow in these basins, the relationship of zero flow days with other hydrologic metrics was explored. Though it may seem counter intuitive to directly compare days with no flow with other metrics that include non-zero flows, these both can be seen as different ways of characterizing the flow regime at a given site. For CV calculations, it was suspected that zero flow days were inflating CV values, and thus also evaluated CV for flow records after removal of zero flow days. POT or R-B Flashiness Index on the non-zero-flow records were not evaluated as POT would not change and R-B Flashiness is calculated by evaluating differences between adjacent timepoints and thus would not be valid in the non-zero records with many missing timepoints.

R (R Core Team, 2018) was used to calculate all hydrologic metrics not generated by the Indicators of Hydrologic Alteration software and for all statistics evaluating potential drivers of observed trends in hydrologic metrics. Prior to statistical analysis, any data that were not normally distributed were transformed using lognormal transformations. Linear mixed effects models (lmer package in R) were used to evaluate relationships among the hydrologic metrics and watershed characteristics, including per cent imperviousness, area-weighted flowrate (q), annual precipitation, watershed area, and watershed slope; both individual variables as well as interactions were included. The models used data from all 14 years of record, where one value was included for each year (sum, mean, or median as described above and was appropriate for the metric). Values were included from each year in order to control for sources of variability during the time period (e.g. precipitation changes year-to-year) but trends over time were not assessed. The study site/watershed was included as a random effect. The Akaike Information Criterion (AIC) was used to compare model results for a given hydrologic metric; the best model was chosen based on the lowest relative AIC. Two watersheds were outliers with respect to the zero flow days metric due to their combination of high average slope, large watershed area, and high precipitation that drove much more consistent flow and lower zero flow days; these watersheds were removed from the zero-flow days linear mixed effects model. Given limited information on stormwater control measures installation year to estimate the per cent stormwater control measures area for each year, a simpler linear regression model was used where each input metric was based on the mean of the 14 year period. To enable comparison between results of this study and those from another geographic region (Mogollon et al., 2016), the relationship between R-B Flashiness and development, rather than imperviousness, was examined.

Baker, D. B., Richards, R. P., Loftus, T. T., and Kramer, J. W. (2004). A new flashiness index: Characteristics and applications to midwestern rivers and streams. Journal of the American Water Resources Association, 40(2), 503–522. https://doi.org/10.1111/j.1752-1688.2004.tb01046.x

Bevelhimer, M. S., McManamay, R. A., and O’Connor, B. (2015). Characterizing Sub-Daily Flow Regimes: Implications of Hydrologic Resolution on Ecohydrology Studies. River Research and Applications, 31(7), 867–879. https://doi.org/10.1002/rra.2781

Homer, C. G., Dewitz, J. A., Yang, L., Jin, S., Danielson, P., Xian, G., et al. (2015). Completion of the 2011 National Land Cover Database for the conterminous United States-Representing a decade of land cover change information. Photogrammetric Engineering and Remote Sensing, 81(5), 345–354.

Mallakpour, I., and Villarini, G. (2016). Investigating the relationship between the frequency of flooding over the central United States and large-scale climate. Advances in Water Resources, 92, 159–171. https://doi.org/10.1016/j.advwatres.2016.04.008

McDaniel, R. D., and O’Donnell, F. C. (2019). Assessment of Hydrologic Alteration Metrics for Detecting Urbanization Impacts. Water, 11(5), 1017. https://doi.org/10.3390/w11051017

Mediero, L., Santillán, D., Garrote, L., and Granados, A. (2014). Detection and attribution of trends in magnitude, frequency and timing of floods in Spain. Journal of Hydrology, 517, 1072–1088. https://doi.org/10.1016/j.jhydrol.2014.06.040

Mogollon, B., Frimpong, E. A., Hoegh, A. B., and Angermeier, P. L. (2016). Recent Changes in Stream Flashiness and Flooding, and Effects of Flood Management in North Carolina and Virginia. Journal of the American Water Resources Association, 52(3), 561–577. https://doi.org/10.1111/1752-1688.12408

Moltz, H. L. N., Palmer, J. B., and Smith, Z. M. (2018). Streamflow Alteration from Impervious Cover: Are All Watersheds Created Equal? JAWRA Journal of the American Water Resources Association, 0(0). https://doi.org/10.1111/1752-1688.12681

R Core Team. (2018). R: A language and environment for statistical computing (Version 3.4). Vienna, Austria: R Foundation for Statistical Computing. Retrieved from https://www.r-project.org/

Schoonover, J. E., Lockaby, B. G., and Helms, B. S. (2006). Impacts of land cover on stream hydrology in the west Georgia piedmont, USA. Journal of Environmental Quality, 35(6), 2123–2131. https://doi.org/10.2134/jeq2006.0113

Smith, B. K., and Smith, J. A. (2015). The Flashiest Watersheds in the Contiguous United States. Journal of Hydrometeorology, 16(6), 2365–2381. https://doi.org/10.1175/JHM-D-14-0217.1

The Nature Conservancy. (2006). IHA Applications. Retrieved from https://www.conservationgateway.org/ConservationPractices/Freshwater/EnvironmentalFlows/MethodsandTools/IndicatorsofHydrologicAlteration/Documents/iha-apps.pdf

The Nature Conservancy. (2009). Indicators of Hydrologic Alteration: Software for Understanding Hydrologic Changes in Ecologically-Relevant Terms. Retrieved February 18, 2019, from https://www.conservationgateway.org/ConservationPractices/Freshwater/EnvironmentalFlows/MethodsandTools/IndicatorsofHydrologicAlteration/Pages/IHA-Software-Download.aspx

US Geological Survey. (2017). USGS-R/EflowStats. R, USGS-R. Retrieved from https://github.com/USGS-R/EflowStats

Xian, G., Homer, C. G., Dewitz, J. A., Fry, J., Hossain, N., and Wickham, J. D. (2011). The change of impervious surface area between 2001 and 2006 in the conterminous United States. Photogrammetric Engineering and Remote Sensing, 77(8), 758–762.