Purpose and Scope

The purpose of this report is to present methods and results for estimating the magnitude and frequency of floods for unregulated streams on five of the main Hawaiian Islands: Kauaʻi, Oʻahu, Molokaʻi, Maui, and Hawaiʻi. (Unregulated streams are those for which peak flows are not altered to a large extent by upstream reservoirs, dams, diversions, or other structures.) The report describes (1) updated flood flows at selected streamgages for the 50-, 20-, 10-, 4-, 2-, 1-, 0.5-, and 0.2-percent AEPs using peak-flow data through water year 2020, (2) regional-skew estimates for Hawaiʻi, (3) regression equations for relating basin and climatic characteristics to flood flows at ungaged stream sites for the 50-, 20-, 10-, 4-, 2-, 1-, 0.5-, and 0.2-percent AEPs, and (4) the accuracy and limitations of the regression equations. The regression equations developed were incorporated in the web-based USGS StreamStats application for Hawaiʻi (Rosa and Oki, 2010). This report supersedes previous flood-frequency reports in Hawaiʻi (the most recent being Oki and others, 2010) by including more recent peak-flow data through water year 2020 (table 2) and by implementing advances in statistical techniques developed after previous reports were published.

Previous Studies

The USGS began collecting streamflow data in Hawaiʻi in 1909 (Yamanaga, 1972, p. 1). The number of active streamgages increased until the late 1960s and decreased steadily since that time. In water year 2020, the USGS had 137 active streamgages in Hawaiʻi (the number of active streamgages includes streamgages that are not included in the present study).

Previous flood studies in Hawaiʻi include descriptive and quantitative investigations related to storm-drainage standards, peak-flow statistics, and (or) regional regression equations at spatial scales ranging from individual floods to statewide studies. For a detailed list and descriptions of past studies related to floods in Hawaiʻi, see Oki and others (2010, app. A). Prior to the study described in this report, the most recent estimates of peak-flow statistics for streams in Hawaiʻi were published in 2010 (Oki and others, 2010), using data from streamgages on Kauaʻi, Oʻahu, Molokaʻi, Maui, and Hawaiʻi through water year 2008. Oki and others (2010) applied the methods described in Bulletin 17B (Interagency Advisory Committee on Water Data, 1982) to estimate AEP flows at 235 streamgages and developed regression equations to relate flood flows to basin characteristics at ungaged locations. Peak-flow estimates from Oki and others (2010) are compared to those generated by this study in the “Comparison of Results with Previous Studies” section.

Description of Study Area

The eight main islands in the State of Hawaiʻi, located in the north Pacific Ocean, trend northwest to southeast and have a total land area of 6,420 square miles (mi²). The majority of each island was formed by shield-stage volcanic eruptions, and the southeastern islands are geologically the youngest. The area considered for flood-frequency analysis in this study (hereinafter referred to as “study area”) includes the five islands of Kauaʻi, Oʻahu, Molokaʻi, Maui, and Hawaiʻi, which have a total land area of 5,904 mi². The maximum altitudes for the islands in the study area are 5,246 feet (ft) for Kauaʻi; 4,046 ft for Oʻahu; 4,948 ft for Molokaʻi; 10,016 ft for Maui; and 13,781 ft for Hawaiʻi (U.S. Geological Survey, 2013). The individual islands in the study area are commonly divided into two physiographic zones, windward and leeward, based on their exposure to the predominant northeasterly winds.

Drainage basins in the study area are characterized by relatively small sizes, amphitheater-shaped valley heads, steep walls, and gently sloping floors (Wong, 1994). In geologically older areas (for example, northern Kauaʻi) with abundant rainfall, erosion and mass wasting have created large valleys and well-defined stream channels; in geologically younger areas (for example, southern Island of Hawaiʻi), valleys are smaller and well-defined stream channels are uncommon because high-permeability soils and rocks at the surface allow rainwater to infiltrate before eroding the land. The topography of Hawaiian shield volcanoes leads to a radial drainage pattern where streams tend to flow away from each other (Oki and others, 2010).

Streamgage Selection and Peak-Flow Data

Peak-flow data for streamgages with at least 10 years of record were downloaded from the USGS National Water Information System (NWIS; https://waterdata.usgs.gov/nwis) database (U.S. Geological Survey, 2021); flood-frequency analyses at a streamgage with less than 10 years of record are generally unreliable (England and others, 2019, p. 36). Data through water year 2020 were used at all streamgages where available. At the time of analysis, USGS streamgages 16103000 (Hanalei River near Hanalei, Kauaʻi) and 16325000 (Kamananui Stream at Pupukea Mil Road Oʻahu) did not have approved 2020 peaks. Additionally, 2021 peaks for USGS streamgages 16345000 (Opaeula Stream near Wahiawa, Oʻahu) and 16587000 (Honopou Stream near Huelo, Maui) were included in the analysis because they were the largest floods on record for each streamgage. Streamgages in Hawaiʻi typically are either continuous-record gages or crest-stage gages. Continuous-record gages record the stage (height) of streamflow at short intervals (for example, every 15 minutes, although the recording interval may be automatically decreased during times of rapidly changing flow at higher stages), whereas crest-stage gages record only the maximum stage of floods above a certain threshold. Discharge is computed from stage measurements at streamgages using a site-specific stage-discharge relation or indirect-measurement methods.

Initially, 268 active and discontinued streamgages were identified as potential streamgages to include in the flood-frequency analysis. The number of usable streamgages was reduced to 238 after reviewing the data (figs. 2–6; table 1.1, app. 1). Streamgages with fewer than 10 years of usable peak-flow data were excluded (including some streamgages where peaks were removed in a low-outlier screening that will be described in section, “Low Outliers Identified with the Multiple Grubbs-Beck Test”). Five streamgages with regulated or diverted flow (NWIS qualification code 6) were excluded (table 1.2, app. 1). Streamgages where discharge was affected to an unknown degree by regulation or diversion (NWIS qualification code 5) were retained except for USGS streamgage 16210500 (Kaukonahua Stream at Waialua, Oʻahu), which was excluded because peak flows may have been substantially affected by regulation. Streamgages with drainage basins that had impervious surfaces covering more than 20 percent of the land were excluded because of the potential effects of impervious surfaces and urbanization on peak flows (see section, “Land Cover”). (The exclusion criterion for urbanized streamgages in the current study—greater than 20 percent impervious surface in the drainage basin—is about equal to the exclusion criterion used by Oki and others (2010) for the 238 streamgages in the current study: greater than 20 percent combined medium- and high-intensity development in the drainage basin.) Additionally, data from 18 streamgages were excluded because of possibly inaccurate rating curves or other potential issues (R.A. Fontaine, U.S. Geological Survey, written commun., 2020). Streamgages and data excluded from the analysis are available in appendix 1 (table 1.2).

From the 238 streamgages used in the study, the average number of available annual peaks for each streamgage is 40. The longest available record, 109 annual peaks, is from Honopou Stream on Maui (USGS streamgage 16587000). The number of available annual peaks for this study reached a maximum of 169 peaks in 1967–68 and has steadily decreased since (fig. 7). About 81 percent of annual peak used in this study with known dates occurred during the wet season from October to April (fig. 8).

This study includes 18 streamgages that were not used by Oki and others (2010)—the most recent flood-frequency study for Hawaiʻi—and excludes 15 streamgages that were used by Oki and others (2010) (tables 1.2 and 1.3, app. 1). Most of the newly included streamgages had fewer than 10 usable annual peaks in 2010; most of the newly excluded streamgages were omitted because they had fewer than 10 usable annual peaks after screening for low outliers.

Trends in Peak Flows

The Bulletin 17C methods for flood-frequency analysis (England and others, 2019) used in this study assume that (1) the peak-flow data are a random, independent, and identically distributed sample that is representative of the population of floods, and (2) the parameters describing the statistical distribution of floods will not change in the future (that is, the distribution is stationary). These assumptions may be violated by deterministic trends related to abrupt or gradual changes in stream regulation, land cover and land use, or climate, or a mixture of those sources (Milly and others, 2008; Vogel and others, 2011). Stationarity can be difficult to detect in hydrologic time series, however, because natural processes often exhibit low-frequency deviations that persist for decades or centuries (Cohn and Lins, 2005; Villarini and others, 2009; Lins and Cohn, 2011). To evaluate the assumptions of the Bulletin 17C methods, peak-flow data were tested for monotonic trends and step trends. Trends were considered statistically significant for probability values (p-values) less than or equal to 0.05—this is the probability that an observed trend is due to random chance.

Peak-Flow Trend Results

Statistically significant monotonic trends were detected at 51 of the 238 streamgages used in this study for at least one of the Mann-Kendall trend-test versions (app. 2). Under the independence assumption, 29 streamgages had decreasing trends and 18 streamgages had increasing trends; under the STP assumption, 26 streamgages had decreasing trends and 16 streamgages had increasing trends; and under the LTP assumption, 6 streamgages had decreasing trends and 4 streamgages had increasing trends. Of the 90 active streamgages (those with annual peaks reported for water year 2020), 24 have significant trends: 16 decreasing and 8 increasing. Of the streamgages with significant monotonic trends, streamgages on Kauaʻi, Oʻahu, and the Island of Hawaiʻi had predominantly decreasing trends, whereas streamgages on Maui and Molokaʻi had mixed results. The magnitude of the trend during the record period at each streamgage was estimated using the Theil slope (also known as the Kendall-Theil robust line) (Helsel and others, 2020, p. 332) (app. 2).

The results from trend testing should be reviewed in the context of spatial and temporal data availability. Test results depend considerably upon the time period analyzed: the available peak-flow data sometimes represent neither long nor concurrent record periods (fig. 9).

In some cases, short-term trends related to natural climate variability can be superimposed on long-term deterministic trends. A flag plot was created by applying the Mann-Kendall trend test (with the independence assumption) to all possible pairs of starting and ending years (minimum 10 annual peak discharges) for each streamgage to determine the percentage of streamgages with statistically significant increasing and decreasing trends for various periods (fig. 10). Record periods ending before about 1980 had predominantly increasing trends; record periods ending after about 1980 had predominantly decreasing trends. Additionally, decadal trends in peak flows were analyzed by applying the Mann-Kendall trend test (with the independence assumption) to each record subperiod spanning 10 years (for example, 1975–84) and containing at least eight annual peaks. The results are shown in figure 11, where the values represent percentages of streamgages with increasing or decreasing trends for the decade, plotted on the mid-decade year (for example, for 1975–84, the mid-decade year would be 1979). Increasing peak-flow trends were common during the decades centered around 1946–48, 1963–64, 1987–88, 2001–03, and 2013–14. Decreasing peak-flow trends were common during the decades centered around 1968–73 and 1993–96.

Significant step trends in peak-flow data were found at 44 streamgages (fig. 12; app. 2). Twenty-six streamgages with significant step trends had peak-flow magnitudes that decreased after the change point; 18 streamgages had peak-flow magnitudes that increased after the change point. Decreasing step trends most commonly have change points clustered from the 1970s to the early 2000s, whereas increasing step trends have a broad distribution of change points spanning the 1930s to 2000s. As with monotonic trends, comparisons of step trends and change points between streamgages should be made with caution because of sometimes inconsistent record periods available for analysis (fig. 9). One of the clearest patterns is that, of the 20 streamgages with statistically significant step trends on Oʻahu, 11 have change points during 1968–76.

Time-series data can possess both monotonic trends and step trends, or one type of trend may mask the other. Thirty-four streamgages had a significant monotonic trend (using the independence assumption) and a significant step trend. Further evaluations of trends in the peak-flow data may consider analyzing the time periods before and after significant change points separately and (or) applying adjustments to the data to account for autocorrelation before using the Pettitt test for step trends.

In summary, trend analyses of the peak-flow data used in this study suggest that decreasing trends are more common than increasing trends. For monotonic trends, the number of significant trends decreases when STP and LTP are accounted for. Monotonic and step trend tests suggest that peak-flow magnitude has generally decreased since about 1980, although the number of streamgages with significant trends may be overestimated by the presence of natural hydroclimatologic fluctuations related to LTP.

Stationarity should be the default assumption for flood-frequency analyses, unless the nonstationarity assumption can be justified based on a clear understanding of the physical processes of trends (Lins and Cohn, 2011; Serinaldi and Kilsby, 2015; Salas and others, 2018; England and others, 2019; Ryberg and others, 2020). Although this study presents a cursory summary of peak-flow trends at streamgages, an exhaustive investigation of trends and the potential causes of trends is beyond the scope of this report. If additional data and more comprehensive analyses discover relations not discussed here, future flood-frequency analyses may consider incorporating trends and nonstationarities and (or) excluding stations with definitive nonstationarities related to deterministic trends. In the absence of clear relations between trends and hydroclimatological forcings, the assumptions in Bulletin 17C (England and others, 2019) are presumed to be valid and are retained for this study.

Physical and Climatic Basin Characteristics

Drainage basins for each streamgage used in the study were delineated using geographic information system (GIS) methods. Physical and climatic basin characteristics for each streamgage were estimated using available datasets. Accurate basin delineations and basin characteristics are critical for the regional regression analysis relating basin characteristics to flood discharges. The source data used to determine the basin delineations and basin characteristics were incorporated into the USGS StreamStats application for Hawaiʻi (Rosa and Oki, 2010). StreamStats allows users to select any point along a stream and automatically delineate a drainage basin and compute selected basin and climatic characteristics for the watershed upstream from that point (U.S. Geological Survey, 2019).

Regional Skew Coefficient

This section presents a general overview of the regional skew coefficient and criteria used for selection of sites in the regional skew analysis. More details about the regional skew regression analysis, including the methodology and calculations, are located in appendix 3.

To help improve estimates of annual peak discharges corresponding to various AEPs—particularly for streamgages with short annual peak-flow records (that is, streamgages with fewer than about 30 annual peaks)—current guidance for flood-frequency analysis by Federal agencies (Bulletin 17C; England and others, 2019) recommends using a weighted average of the at-site and regional skews. Previous guidance (Bulletin 17B; Interagency Advisory Committee on Water Data, 1982) supplied a national map of regional skew but encouraged hydrologists to develop more localized models when appropriate.

Because of complications introduced using EMA and MGBT (Cohn and others, 1997) and large cross-correlations between annual peak discharges at pairs of streamgages, a Bayesian weighted least-squares/Bayesian generalized least-squares (B–WLS/B–GLS) regression framework was developed to provide stable and defensible results for regional skew (Veilleux, 2011; Veilleux and others, 2011; Lamontagne and others, 2012; Veilleux and others, 2012). B–WLS/B–GLS uses ordinary least-squares (OLS) regression to fit an initial model of regional skew that is used to generate a stable estimate of regional skew for each streamgage. This estimate is the basis for computing the variance of each estimate of at-site skew used in the B–WLS analysis. B–WLS is then used to generate estimators of the regional skew model parameters. Finally, B–GLS is used to estimate the precision of those estimators, the model error variance and its precision, and compute various diagnostic statistics.

In this study, EMA with MGBT was used to estimate the at-site skew, G, and its mean squared error, MSE_G. EMA with MGBT allows for the censoring of low floods as well as the use of flow intervals to describe missing, censored, and historical data. EMA with MGBT complicates the calculations of effective record length (and effective concurrent record length) used to describe the precision of skew estimates because the annual peak discharges are no longer represented by single values. To properly account for these complications, the B–WLS/B–GLS procedure was used.

A total of 124 streamgages that were not redundant (see section, “Elimination of Redundant Sites”) and had a pseudo effective record length (P_RL) of 36 years or greater were used to develop the final regional skew model for Hawaiʻi (table 3.1, app. 3; for more information on pseudo effective record length, see app. 3). To explain the variability in skew, a windward/leeward split of flood regions and 17 basin characteristics were tested, but this approach failed to provide sufficient predictive power. Therefore, a constant regional skew (G_R) of −0.157 was selected for the State of Hawaiʻi (app. 3). The average variance of prediction (AVP_new), 0.212, is equivalent to the mean square error of the regional skew (MSE_R) and corresponds to an effective record length of 36 years. These values supersede G_R (−0.05), MSE_R (0.302), and effective record length of 17 years associated with the generalized skew map in Bulletin 17B (Interagency Advisory Committee on Water Data, 1982), which was used in the previous study (Oki and others, 2010).

Because of the relatively large uncertainty in the at-site skew for short to modest record lengths, the at-site skew and its MSE can be weighted with the regional skew and its MSE to generate a better, weighted estimate of skew for a given streamgage basin (Tasker, 1978; England and others, 2019, app. 7). Large deviations between the at-site and regional skew may indicate that the flood frequency characteristics of the basin of the streamgage of interest differ from those used to estimate the regional skew. If the at-site and regional skews differ by more than 0.5, it is considered reasonable to use the at-site skew instead of the weighted skew in the EMA (England and others, 2019, p. 25–26). The weighted skew was used within PeakFQ at all except the following streamgages, where the at-site skew was used: USGS streamgages 16400000 (Halawa Stream near Halawa, Molokaʻi), 16501200 (Oheo Gulch at dam near Kipahulu, Maui), 16502000 (Hahalawe Gulch near Kipahulu, Maui), 16557000 (Alo Stream near Huelo, Maui), 16565000 (Kaaiea Gulch near Huelo, Maui), 16638500 (Kahoma Stream at Lahaina, Maui), and 16717600 (Alia Stream near Hilo, Island of Hawaiʻi).

Expected Moments Algorithm Frequency Analysis

The guidelines in Bulletin 17C (England and others, 2019) suggest using the expected moments algorithm (EMA) to analyze the available flood data. EMA improves upon the methods provided in the previous flood-frequency guidelines, Bulletin 17B (Interagency Advisory Committee on Water Data, 1982), by cohesively incorporating all available flood-related data, including historical flood information, zero flows, low outliers, flow intervals, and perception thresholds (Lane and Cohn, 1996; Cohn and others, 1997; England and others, 2019).

Flood-frequency data generally come from two types of sources: systematic and historical data. Systematic data are the primary source of flood-frequency data for Hawaiʻi and consist of peak-flow data collected at regular intervals from either continuous-record gages or crest-stage gages. Systematic data are usually collected in consecutive years, although records sometimes contain data gaps between years (for example, see fig. 9). Historical data consist of major floods that exceeded a perception threshold and occurred outside the period of routine streamgaging, independent of how recently the flood occurred. Historical floods are valuable because they can be used to extend records with the knowledge that if a particular discharge was exceeded, it would have been recorded in some way. For example, at Kalihi Stream (USGS streamgage 16229300), historical data indicate that the flood on May 14, 1960, was the largest flood since at least 1937; with this information, the EMA can incorporate the period 1937–59 into the analysis by indicating that all annual peaks during this period were less than the peak discharge on May 14, 1960: 6,350 ft³/s.

Some peaks in the peak-flow record are classified as “opportunistic.” Opportunistic peaks occurred outside the period of systematic streamgaging and were measured because of operational decisions other than the exceedance of a perception threshold. Because the statistical sampling properties of opportunistic peaks are unknown, opportunistic peaks were excluded from flood-frequency analysis.

Low Outliers Identified with the Multiple Grubbs-Beck Test

Peak-flow records commonly contain low-magnitude outliers that deviate considerably from the rest of the peak population. Low outliers often have a disproportionately large influence on the fit of the frequency curve, at the expense of the fit at the high-discharge end of the curve (Cohn and others, 2013). Because most applications of flood-frequency analysis (for example, infrastructure design and flood protection) focus on the lower AEPs (larger peak flows), low outliers are removed when fitting the frequency curve. Bulletin 17C guidelines (England and others, 2019) recommend the Multiple Grubbs-Beck Test (MGBT) to detect and remove the potentially influential low floods (PILF) during flood-frequency analysis. The MGBT improves upon the Grubbs-Beck test (Grubbs and Beck, 1972) recommended in Bulletin 17B by accommodating the possibility that several low floods are potentially influential (Cohn and others, 2013; Lamontagne and others, 2016). The PeakFQ program, version 7.3, automatically applies the MGBT when computing flood statistics. For a few streamgages, a threshold for PILF detection and removal was manually set to improve the fit at the upper end of the frequency curve. For an example of a fitted frequency curve where PILFs have been removed, see figure 14.

Flood-Frequency Estimates at Gaged Sites

Flood-frequency estimates for 238 streamgages in Hawaiʻi were calculated using EMA and MGBT techniques and the new regional skew coefficient (see section, “Regional Skew Coefficient”). The magnitude of peak flows for the 50-, 20-, 10-, 4-, 2-, 1-, 0.5-, and 0.2-percent AEPs are listed in appendix 4. Three estimates of peak flows for the select AEPs are provided: the at-site estimate (EMA), the regression-equation estimate (regression), and a weighted average of the other two estimates (weighted). The weighted estimate is generally preferred for most situations because it combines information from the independent at-station and regression-equation estimates (England and others, 2019, p. 33). The regression and weighted estimates will be discussed in subsequent sections.

Elimination of Redundant Sites

Streamgages on the same stream and with similar size drainage basins may contain redundant hydrologic information. Redundant streamgages generally have similar basin characteristics and hydrologic responses to a given storm; thus, they provide only one unique statistical observation, rather than two independent observations. Including redundant streamgages in a regression analysis can negatively affect the results because the data are not independent (Farmer and others, 2019). To determine if two streamgages should be classified as redundant, three types of information are reviewed: (1) whether the streamgages have nested drainage basins, meaning that one drainage basin is entirely contained within the other, (2) whether the streamgages have drainage basins of similar size, and (3) whether the streamgages have temporal overlap in their peak-flow data.

To evaluate the likelihood that the streamgages have nested drainage basins, the standardized distance was computed. The standardized distance (SD) between two basin centroids is defined as:

A drainage area ratio was used to determine if the two streamgages had similarly sized drainage areas. The drainage area ratio (DAR) is defined as:

A script written in the R programming language (R Core Team, 2021) was used to provide an initial screening of potentially redundant streamgages on the basis of the standardized distance and drainage area ratio. Screening thresholds for the standardized distance and drainage area ratio were set to 0.5 and 5, respectively. All possible combinations of streamgage pairs from the 238 streamgages were considered in the redundancy analysis. The script identified 50 potentially redundant pairs of streamgages (some streamgages were identified as potentially redundant to more than one other streamgage); 9 of the streamgage pairs identified were not nested or did not have temporal overlaps in their data and therefore were not considered redundant. For streamgage pairs that had nested basins and drainage area ratios less than five, one of the streamgages was classified as redundant and removed from the regression analysis. Generally, streamgages with longer periods of record were prioritized and retained. Additional considerations include (1) whether the streamgage is active and (2) how the streamgage’s basin characteristics fit into the distribution of basin characteristics for the streamgages in the region. For example, the streamgage with the largest drainage basin in the region may be prioritized, all other factors being equal, because it expands the range of values used to develop the regression equations (as will be discussed later, the regression equations should not be used with values for the explanatory variables that are not within the range of values used to develop the regression equations). Of the 238 streamgages considered for the regression analysis, 23 streamgages (about 10 percent) were classified as redundant and removed from the analysis (table 1.1, app. 1). Redundant streamgages were not included in either the regression or skew analyses, although station-specific flood statistics for redundant streamgages were still computed (table 1.1, app. 1).

Exploratory Data Analysis

An exploratory data analysis was completed to evaluate the best combinations of basin characteristics to use as explanatory variables in the regional regression equations. Scatter plots of each pair of response variable (flood discharge associated with an AEP) and explanatory variable (basin characteristic) were created to visually evaluate the pattern and linear relation between variables. Because multiple-linear regression seeks to quantify linear relations between the response and explanatory variables, data transformations were tested to improve the linearity of the relation. Logarithmic transformations (base 10) of both response and explanatory variables resulted in the best linear relations and the most constant variance (homoscedasticity) about the regression line. Prior to log transformation, a constant of 1 was added to the explanatory variables expressed as a percentage (for example, percentage impervious land cover) because some values for these basin characteristics are 0, and the logarithm of 0 is undefined. When logarithmic transformations are applied to the data to improve linearity, equation 2, which describes the general form of the regression model, becomes:

OLS regression analysis was performed using the R programming language (R Core Team, 2021) and the “smwrStats” statistical package (U.S. Geological Survey, 2017). An “all-possible subsets” regression was used to determine the best combinations of explanatory variables for models with 1–3 variables. Generally, the regression equations were limited to 1 variable per 10 streamgages in the region (Farmer and others, 2019). Model diagnostics for regression equations were reviewed for adequacy; selection criteria for explanatory variables included: (1) maximize the adjusted coefficient of determination, (2) minimize Mallow’s C_p statistic, and (3) minimize the predicted residual sum of squares (PRESS) statistic. Basin characteristics and the magnitude and sign of their coefficients in the regression equations were reviewed to ensure hydrologic plausibility. Multicollinearity, where two explanatory variables have a strong linear dependency, was evaluated using the variance inflation factor (VIF)—no variables in the final OLS regression models had a VIF greater than 5, indicating that multicollinearity is unlikely.

Regional Regression Equations

GLS multiple-linear regression (Stedinger and Tasker, 1985; Tasker and Stedinger, 1989; Farmer and others, 2019) was used to determine the final coefficients and performance metrics for each regional regression equation. GLS regression techniques are generally preferred for flood-frequency analysis because they improve estimates of AEP discharges and estimates of the accuracy of the regression model by accounting for (1) unequal record lengths from streamgages and (2) cross-correlation of streamflow statistics from streamgages (England and others, 2019, p. 30–31). Streamgages with shorter records are given less weight than streamgages with longer records. Additionally, less weight is given to streamgages where concurrent peak flows are correlated with nearby streamgages. Based on results from the OLS analysis, all explanatory and response variables were log-transformed prior to use in the GLS regression analysis. The USGS weighted-multiple-linear regression (WREG) program version 2.02 (https://github.com/USGS-R/WREG)—written in the R programming language (R Core Team, 2021)—was used to compute the final GLS equations and model performance metrics. Default values of the correlation-smoothing function in WREG (alpha, α, 0.002; and theta, θ, 0.98) were adjusted for each regional regression equation to improve model fit.

The final regression equations for estimating AEP discharges from basin characteristics for ungaged streams in Hawaiʻi are shown in table 5. Within each region, the explanatory variables were constrained to be identical for all regression equations; this ensures that predicted discharges uniformly increase as the AEP decreases. All explanatory variables used in the regression equations were statistically significant at the 95-percent confidence level (p-value less than or equal to 0.05) for at least one AEP in each region. Ranges of the explanatory variables used to develop the regression equations are presented in table 6.

Drainage area is the most common explanatory variable and appears in all equations except for those representing southern Island of Hawaiʻi (region 10). Other explanatory variables used in the regression equations include mean annual precipitation, precipitation-frequency statistics (for example, the maximum 48-hour rainfall that occurs on average once in 500 years, or I48H500Y), and soil permeability. The only region containing more than 1 explanatory variable per 10 streamgages was leeward Kauaʻi (region 1); here, a second explanatory variable—mean annual precipitation—was justified for the 17 streamgages in the region because of its clear hydrologic relation to flood magnitude and the notable improvement in performance metrics of the regression model, relative to a one-variable (drainage area) equation. All explanatory variables used in the regression equations except for soil permeability had positive coefficients, indicating that as the value of explanatory variable increases, the predicted flood discharge also increases.

Soil permeability was used as an explanatory variable in the leeward Oʻahu regression equations (region 3). The negative coefficient of the variable in all equations indicates that, as soil permeability increases, the predicted peak discharge decreases; this is intuitive because more permeable soils allow more water to infiltrate and reduce the volume of overland runoff than can contribute to flood peaks. The statistical significance of soil permeability decreases at the low-AEP floods and the p-value for soil permeability as a variable exceeds 0.05 at an AEP of 0.002. The importance of soil permeability is expected to be lower for the low-AEP floods because most soils become fully saturated during large flooding events; thus, the spatial variability of permeability will have a smaller influence on the flood magnitude above a certain threshold (Hollis, 1975; Konrad, 2003).

Southern Island of Hawaiʻi (region 10) is the only region where drainage area was not used as an explanatory variable. This may reflect the lack of well-defined stream channels in this area and the high permeabilities of soils and rocks at the surface (Oki and others, 2010). Because rainfall tends to infiltrate the surface quickly, drainage area likely does not play a large role in governing flood magnitudes for most AEPs in region 10.

Leverage and influence statistics computed by WREG were reviewed to evaluate how each streamgage affected the regression results. Leverage is a measure of how far away the values from one streamgage are from the values at all other streamgages in the regression model and is used to identify unusual observations (Farmer and others, 2019). Streamgages with high leverage have the potential to exert a strong influence on the regression parameters (Helsel and others, 2020, p. 238). The influence metric is a measure of how much influence a particular streamgage has on the regression parameters. Streamgages with leverage and influence values that exceeded the thresholds calculated by WREG were reviewed for potential errors in the peak-flow data (for example, poor computation of discharge from stage or inaccurate data representation in NWIS) and basin-characteristic data (for example, inaccurate basin delineations) or other issues that would make the streamgage ineligible for the regression model. Although several streamgages had both high leverage and influence, no errors were identified in the associated data and a reasonable hydrologic justification for removing the streamgages from the regression analysis could not be found.

Assessment of Fit

To assess the fit of the regression models, peak discharges predicted using the regression equations can be compared to the peak discharges determined from fitting the observed peak-flow data to the LP-III frequency curve (fig. 15). Comparisons of the observed and predicted peak discharges suggest that the regression models fit the observed data well. Additionally, plots of residuals (the difference between predicted and observed discharges) were examined to test the validity of the assumptions related to the regression models (fig. 16). The plots show that the residuals do not follow any trend and are equally distributed around zero and indicate that no model assumptions were violated.

Accuracy and Limitations of Regional Regression Equations

Several performance metrics from the WREG program can be used to evaluate the accuracy of the regression equations. Performance metrics for the final regression equations are presented in table 5 and include the model error variance (MEV), the standard model error variance (SMEV), the root mean squared error (RMSE), the pseudo coefficient of determination (pseudo-R²), the average variance of prediction (AVP), and the average standard error of prediction (SEP_avg). The MEV and SMEV describe the portion of the total error that can be attributed to having an imperfect model, in log units and percent, respectively. The RMSE and pseudo-R² are measures of model accuracy for streamgages used in the model development. RMSE describes how much the predicted peak discharges deviate from the observed peak discharges. The pseudo-R² is a measure of the predictive strength of the regression model and describes the variability of the response variable that is explained by the explanatory variables, after accounting for the effect of time-sampling error. The pseudo-R² is similar to the standard coefficient of determination (R²), where the closer the value is to 1.0 (or 100 percent), the greater the amount of variance that is explained by the regression. The AVP and SEP_avg are measures of how well the regression model performs at predicting peak discharges for ungaged sites not used to develop the regression equations—lower values indicate greater predictive power. Generally, pseudo-R² and AVP were the most important metrics for selecting the final regression equations. Equations for calculating these performance metrics are available in Eng and others (2009).

The pseudo-R² values for the final regression equations ranged from 11.9 percent (region 7; 0.5 AEP) to 100.0 percent (region 6; 0.002–0.02 AEPs) for all regions and AEPs, and ranged from 46.2 percent (region 5) to 100.0 percent (region 6) for the 0.01-AEP discharges. The AVP values for the final regression equations ranged from 0.0046 (region 6; 0.02 AEP) to 0.3400 (region 5; 0.5 AEP) for all regions and AEPs, and ranged from 0.0062 (region 6) to 0.2455 (region 10) for the 0.01-AEP discharges. Overall, the final regression models performed well for all regions except regions 5 and 10. The large errors associated with the regression equations for regions 5 (leeward Molokaʻi) and 10 (southern Island of Hawaiʻi) reflect the relatively few sites and short records available and the poor understanding of flood-producing mechanisms in these areas. Collection of additional peak-flow data and consideration of basin characteristics that better characterize flood conditions may improve the accuracy of regression equations used to predict peak flows at ungaged sites in these regions.

The accuracy of the regression equations may be affected by issues associated with the explanatory and response variables. The basin characteristics used as explanatory variables rely on accurate basin delineations—any change to the basin delineations will affect the computed values for each basin characteristic, even if the source data for the basin characteristics remain unchanged. Future analyses that use higher-resolution DEMs (for example, DEMs derived from light detection and ranging [lidar] data) to delineate basins will result in improved accuracy for basin characteristics used in the regression analyses. Potential errors associated with the peak-flow statistics used as response variables in the regression equations include non-representative data from streamgages with short records and errors computing discharge from stream stage (Turnipseed and Sauer, 2010).

The explanatory variables (basin characteristics) should be representative of basin conditions during the period of record that was used to derive the response variables (peak-flow statistics) (Farmer and others, 2019). The historical periods used to determine the mean annual precipitation, precipitation-frequency statistics, and soil permeability used in the final regression equations were 1978–2007, 1899–2005, and 2001–19, respectively. The peak-flow statistics used as response variables were determined from peak-flow data collected during 1911–2021; however, the temporal availability of peak-flow data for each region varies substantially (for example, see fig. 9). If a streamgage used in the regression analysis only has peak-flow data collected during a particularly wet period (for example, the 1960s), it may not accurately characterize the expected long-term relation between basin characteristics and peak-flow statistics in the region.

The regression equations are only applicable to ungaged locations with basin characteristics within the range of values used to develop the regression equations (table 6). If the regression equations are used beyond these limits, the accuracy of the estimated peak-flow statistics would be unknown. Additionally, basin characteristics at ungaged sites should be computed using the same datasets and methods as were used in this study. The USGS StreamStats web application provides the same datasets and methods as were used in this study and therefore provides the most consistent way to estimate peak-flow statistics at ungaged sites (see section, “Estimating Flow Statistics Using StreamStats”). The regression equations developed in this study apply to stream sites that are unregulated and have less than 20 percent impervious land cover; however, StreamStats is unable to warn users or prevent them from selecting regulated stream sites because identifying all stream sites in the study area that are substantially affected by regulations is challenging.

Uncertainty of Individual Estimates Computed Using the Regression Equations

When applying the regression equations at ungaged sites not used in the development of the regression equations, it is important to understand the accuracy and uncertainty associated with the estimated discharge. Two commonly used metrics for the accuracy of an estimated discharge at a particular site are the variance of prediction (VP) and the standard error of prediction (SEP). For a site i, the variance of prediction is the sum of the model error variance and the sampling error variance. The variance of prediction for each streamgage used in the development of the regression equations is computed in WREG. To compute the variance of prediction for estimated discharges using the regression equations ($V{P}_{reg}$), the following equation can be applied:

The sampling error variance is computed as follows:

The SEP for site i, expressed in log units, can be computed as:

To understand the uncertainty of an individual discharge estimate when using the regression equations at an ungaged site, prediction intervals can be computed. Prediction intervals combine the uncertainty of the regression parameters and placement of the regression line—typically described by confidence intervals—with the uncertainty associated with the residuals (Farmer and others, 2019). Prediction intervals describe the range of values within a specific confidence interval, within which the true value exists (Helsel and others, 2020). The following equations can be used to compute the 95-percent prediction intervals for peak-flow estimates obtained using the regression equations for ungaged locations:

For sites not used in the development of the regression equations, the average variance of prediction (AVP), available in table 5, should be substituted for the site-specific variance of prediction ($V{P}_{re{g}_i}$) in equations 9 and 10.

Weighting Flood-Frequency Estimates at Gaged Sites

An improved estimate of peak discharge can be obtained for gaged sites by combining information from the at-site flood-frequency curve (that is, EMA) with information from the regional regression equations (England and others, 2019). The two estimates of peak discharge can be considered independent if the regression equations were developed with a large number of sites. The combined peak-flow estimate is weighted based on the inverse of the variance of prediction for each independent estimate. The weighting procedure only applies to unregulated streams with minimal basin urbanization. The weighted discharges are available in appendix 4 (and can be accessed using StreamStats) and were computed with the following equation:

The variance of prediction for the weighted-discharge estimate ($V{P}_{wt{d}_i}$) can be computed as follows:

Once $Q_{wt{d}_i}$ and $V{P}_{wt{d}_i}$ have been determined, the upper and lower prediction intervals for the weighted-discharge estimate can be computed as follows:

Weighting Flood-Frequency Estimates at Ungaged Sites with Data from a Nearby Gage

Peak-flow estimates determined using the regional regression equations for ungaged sites can be improved if the ungaged site is located near a gaged site on the same unregulated stream. This approach combines information from the regression equations at the ungaged site with information from the discharge estimates at a nearby gaged site by, in part, relating the drainage areas of the two sites. The weight of the regression-derived discharge estimate for the ungaged site increases relative to the weight of the discharge estimates from the gaged site as the drainage-area ratio increases (that is, as the distance between the two sites increases). To apply this weighting method, the gaged site must have at least 10 years of peak-flow data, and the ungaged site must have a drainage area within 50 percent of the drainage area of the gaged site (the drainage-area ratio is more than 0.5 and less than 1.5) and must be located on the same stream (Ries, 2007; Feaster and others, 2009). The following equation can be used to compute the weighted-discharge estimate at the ungaged site: