Purpose and Scope

The purpose of this report is to present methods for estimating the magnitude of peak flows and low flows at defined frequencies and durations on nontidal streams in Delaware. This report (1) describes methods used to estimate the magnitudes of both peak flows and low flows at sites monitored by streamflow-gage sites; (2) presents estimates of peak-flow magnitude at the 50-percent (2 year [yr]), 20-percent (5 yr), 10-percent (10 yr), 4-percent (25 yr), 2-percent (50 yr), 1-percent (100 yr), 0.5-percent (200 yr), and 0.2-percent (500 yr) annual-exceedance probability (AEP) for 94 streamflow-gage sites in and near Delaware; (3) presents estimates of low-flow magnitudes at selected frequencies (7Q2, 7Q10, 7Q20, 14Q2, 14Q10, 14Q20, 30Q2, 30Q10, and 30Q20) describing 7-, 14-, and 30-consecutive-day low-flow periods with annual nonrecurrence intervals of 2, 10, and 20 years (annual non-exceedance probabilities of 50 percent, 10 percent, and 5 percent, respectively) for 45 streamflow-gage sites in and near Delaware; (4) describes methods used to develop regression equations to estimate the magnitude of peak flows and low flows for defined frequencies and durations at ungaged sites in Delaware; (5) describes the accuracy and limitations of the peak-flow and low-flow equations; (6) presents example applications of the peak-flow and low-flow methods; and (7) describes the U.S. Geological Survey (USGS) StreamStats web application from which basin characteristics, streamflow statistics, and estimates can be obtained.

Previous Investigations

Methods for determining peak-flow magnitude estimates for nontidal streams in Delaware were published in previous reports by Tice (1968), Cushing and others (1973), Simmons and Carpenter (1978), Dillow (1996), and most recently by Ries and Dillow (2006). The regionalization methods described in the earlier reports relied on fewer monitored sites and shorter periods of recorded data than those used in Ries and Dillow (2006). Since Ries and Dillow (2006) was published, an additional 13 water years of streamflow data and computationally improved regionalization techniques have become available.

Methods for determining low-flow statistics for nontidal streams in Delaware were published previously in Carpenter and Hayes (1996). Regionalization methods used in that report relied on data through March 1987 (the end of the 1986 climatic year). An additional 31 years of streamflow data and computationally improved regionalization techniques for low-flow statistics have become available since the analyses described by Carpenter and Hayes (1996).

Doheny and Banks (2010) also presented selected low-flow statistics for 114 continuous-record streamflow-gage sites in Maryland based on available data through the 2009 climatic year. Low-flow statistics at several streamflow-gage sites from the Coastal Plain area of Maryland and the Piedmont area of northeastern Maryland included in that investigation were updated using an additional 8 years of continuous streamflow data and used for regionalization in this investigation.

Determination of Basin Characteristics

Basin characteristics that were initially considered for use in the regression analyses included data variables related to drainage area, basin slope, impervious area, soil types, basin storage, developed land, mean basin elevation, basin relief, forest cover, mean-annual and maximum precipitation totals, hydraulic conductivity, development density, outlet-point elevation, flowpath lengths, housing units and housing density, and basin population. A list of specific basin characteristics considered for use in both the peak-flow and low-flow regression analyses is shown in table 2.

Flow Statistics for Streamflow-Gage Sites

There were 121 streamflow-gage sites identified either within the State of Delaware or within a 25-mi buffer of the State of Delaware border and located within Maryland, New Jersey, or Pennsylvania, each of which is associated with 10 or more years of annual peak-flow data. Of these sites, 94 were subsequently used for the development of peak-flow regression equations following analyses of trend, redundancy, and regulation as described in the sections below. Following Wagner and others (2016), record lengths from individual gaged sites were not combined following the analysis of contributing-area redundancy and flow regulation. Record combination would only have been possible for two sets of gage sites, and instead of introducing potential error, the longer record was chosen for each set of redundant gages (01467087 and 01467089, 01476480 and 01476500) providing a dataset composed solely of individual streamflow-gage-site data records.

For each of the 121 streamflow-gage sites, peak-flow magnitudes at the 50-, 20-, 10-, 4-, 2-, 1-, 0.5-, and 0.2-percent AEPs (2-year to 500-year recurrence intervals) were determined using the PeakFQ program version 7.3 (Veilleux and others, 2014) and guidance from Bulletin 17C Guidelines for determining flood-flow frequency (England and others, 2019). Bulletin 17C is implemented in PeakFQ version 7.3 using the Expected Moments Algorithm (EMA) (Cohn and others, 1997) and a generalized version of the Grubbs-Beck test for low outliers (Cohn and others, 2013). Flow statistics are provided for sites included in regression analyses as well as those not suitable for inclusion in regression equations owing to regulation, redundancy, or the presence of trend as detailed below. At-site, regression-based, and weighted flow statistics are provided for sites included in regression equation development. At-site, regression-based, and weighted peak-flow estimates along with basin characteristics can be found in Hammond (2021a).

Procedures for Implementing the Bulletin 17C Guidelines

Procedures for implementing the Bulletin 17C guidelines include (1) EMA analysis for fitting the log-Pearson Type III (LP3) distribution, incorporating historical information where applicable; (2) the use of weighted skew coefficients (based on weighting at-site skew coefficients, fig. 3, with regional skew coefficients from the section “Regional Skew Analysis and Determination of Streamflow Analysis Regions” below); and (3) the use of the Multiple Grubbs-Beck Test (MGBT) for identifying potentially influential low flows (PILFs) (England and others, 2019; Sando and McCarthy, 2018).

The standard procedures for incorporating historical information with the EMA methodology involve defining and applying the best-available flow intervals and perception thresholds for peak-flow magnitude analyses that include recorded streamflow data, historical peak flows, and ungaged historical periods. The EMA method improves upon the standard LP3 method in Bulletin 17B by integrating censored and interval peak-flow data into the analysis (Cohn and others, 1997). There are two types of censored peak-flow data in this report: (1) annual peak flows at crest-stage gages for which the flow is less than a minimum recordable value, and (2) historical annual peak flows that have not exceeded a recorded historical peak (Eash and others, 2013). In EMA, interval discharges were used to characterize peak flows known to be greater or less than a specific value, peak flows that could only be reliably estimated within a certain range, or to characterize missing data in periods of recorded data. The MGBT is used by EMA to identify not just low outliers but also other PILFs (Cohn and others, 2013). Although low outliers are typically one or two homogeneous values in a dataset that do not conform to the trend of the other observations, PILFs have a magnitude that is much smaller than the flood quantile of interest, occur below a statistically significant break in the flood magnitude-frequency plot, and have excessive influence on the estimated magnitude of large floods at defined frequencies. An example implementation of Bulletin 17C methods is shown in figure 4 for USGS gage 01480000 (Red Clay Creek at Wooddale, Delaware).

For crest stage gage locations, which differ from continuous record sites in that measurements are only available for the maximum stage of an individual event, it was also necessary to determine the minimum recordable flow as part of implementing perception thresholds. The minimum recordable flow threshold primarily affects treatment of peak flows near the extreme lower tail of the frequency distribution. Generally, the handling of peak flows near the extreme lower tail of the frequency distribution has limited effect on magnitude-frequency analyses, with very low peak flows either censored by PILF thresholds or exerting small influence in determining the distributional parameters of the LP3 distribution.

For some streamflow-gage sites, the peak-flow records are not well represented by the standard procedures and require informed user adjustments. For this study, adjustments were made to account for atypical lower-tail peak-flow data at seven sites (01412500, 01467086, 01467317, 01467351, 01473470, 01475019, and 01488500) by implementing the Single Grubbs-Beck Test instead of the recommended MGBT, resulting in an improved EMA fit of the LP3 distribution. The PeakFQ models for all 121 sites used in this study are provided in Hammond (2021b).

Analysis of Trends in Annual Peak-Flow Time Series

PeakFQ implements the nonparametric Mann-Kendall test for the identification of monotonic trends in annual peak-flow data. The results from this test show that peak-flow data from 18 of the original 121 sites had significant trends, where significance is defined by p-values less than 0.05. Similar to Wagner and others (2016), if sites displayed a significant trend, the Mann-Kendall trend test was implemented on variations of the original annual peak-flow series by removing 5 percent of the annual peak-flow data from the beginning or end of the dataset, or a combination of the two. This practice assumes that if a significant trend is not identified using 90 or 95 percent of the data, then the trend was likely not present and not impactful on the peak-flow magnitude analysis (Eash and others, 2013). If a trend remained after removing data as described, the site was not included in the development of regional regression equations. Nine of the 18 sites were kept for further analysis, although some not included were possibly eliminated owing to redundancy or regulation as described in the sections below. The nine sites that were not used for regression equation development owing to trends were 01467317, 01475000, 01477800, 01484100, 01484270, 01484500, 01491000, 01492500, and 01493000.

Regional Skew Analysis and Determination of Regression Regions

To assess whether separate regions with different regional skew values were present within our study area, we implemented the nonparametric Wilcoxon rank sum test for significant differences between the Coastal Plain and Piedmont physiographic regions used in the prior peak-flow frequency study (Ries and Dillow, 2006). For this analysis, we only considered sites with more than 35 years of peak-flow record, following Bulletin 17B (Interagency Advisory Committee on Water Data, 1982) and Wagner and others (2016). The Wilcoxon rank-sum test showed no significant difference between the long-term skew values in the Coastal Plain and Piedmont physiographic regions; thus, the average of all Coastal Plain and Piedmont at-site skews were used to calculate a regional skew value (0.264) applicable to the entire study region. This value was used along with the at-site skew coefficients to produce the weighted-skew PeakFQ models. Though only one regional skew value was used across the study area, when regression equations were constructed, separate regression equations were built for the Coastal Plain and Piedmont physiographic regions owing to improved model performance for individual regions as compared to using one model for the entire study area. Sites were assigned to physiographic regions based on the region containing the majority of their contributing drainage area. For one site, we deviated from this rule. USGS 01479000 White Clay Creek near Newark, DE, has 55 percent of the contributing area in the Coastal Plain and 45 percent from the Piedmont. The 45 percent located in the Piedmont is in the headwater of the basin where much of the flood generating urban influence is likely located, and this was the justification for assigning it to the Piedmont region in the prior study (Ries and Dillow, 2006).

Sample Adjustments and Removal of Sites Owing to Redundancy

Streamflow-gage sites with significant degrees of regulation were not included in the original dataset if there were regulation qualifiers applied to the annual peak flows in the USGS National Water Information System (NWIS) database. Adjustments to the dataset were also made based on any redundancy among sites that could influence the relations determined by regression analysis. The process used to justify removal of streamflow-gage sites as a result of redundancy is discussed below.

Statistical redundancy occurs when data from one streamflow-gage site can be partly or completely deduced from the data associated with another nearby gage on the same stream or river, or nearby in the same watershed. Inclusion of redundant data in a peak-flow analysis can influence the resulting regression-derived relation. The peak-flow dataset for this investigation was evaluated for redundancy in order to remove any potential influence or bias on the peak-flow magnitude regression analyses. Sites were flagged for possible redundancy if the drainage area ratio between two sites was less than 5.0 and the standardized distance between the two sites was less than 0.5 (Wagner and others, 2016).

There were 25 potential instances of redundancy identified in the peak-flow dataset, involving a total of 30 sites. Decisions regarding site redundancy for peak flows were made based primarily on record length and drainage-area size. Consideration was also given to the length of the data record, any breaks in the continuous record or annual peak-flow record, and how well the peak-flow record represents the current watershed conditions.

Of the 30 sites potentially affected by redundancy issues and questions, 18 were removed from the peak-flow analysis. For this study, no site records for peak flow were combined following removal of regulated or trend-affected sites, as a minimum of 10 years of concurrent data were not available (following Bulletin 17C guidelines and Wagner and others [2016]). The final number of streamflow-gage sites used for peak regression equation development in each state and region is provided in table 3, and a summary of the drainage areas for these streamflow-gage sites is provided in table 4. The ranges of basin characteristics used to develop the equations for each subregion are provided in table 5.

Development of Regression Equations

In developing regional regression equations, basin characteristics for candidate regression models were evaluated based on multicollinearity with other basin characteristics and statistical significance using ordinary least-squares regression. Subsequently, generalized least-squares regression was used to develop the final regression models for each region.

Generalized Least-Squares Regression for Peak Flow Estimation

Generalized least-squares (GLS) multiple-linear regression was used to compute the final regression coefficients and measures of accuracy for the regression equations to estimate peak flows at defined frequencies. The GLS method weights data from streamflow-gage sites in the regression analysis according to differences in streamflow reliability (record length), variability (record variance), and spatial cross correlations of concurrent streamflow among streamflow-gage sites (Farmer and others, 2019). Compared to OLS regression, the GLS regression approach provides improved estimates of streamflow statistics and increases the predictive accuracy of the regression equations (Stedinger and Tasker, 1985). The weighted-multiple-linear regression (WREG) program (Eng and others, 2009) updated to run in R (Farmer, 2017) was used to perform the GLS regressions. A correlation smoothing function is used by WREG to compute a weighting matrix for the streamflow-gage sites in each subregion. The smoothing function relates the correlation between annual peak discharges from each pair of streamflow-gage sites to the geographic distance between the sites. Smoothing functions were developed iteratively and generated separately for the Coastal Plain and Piedmont subregions. The statistical parameters needed to perform prediction-interval computations for peak-flow magnitude estimates can be found in table 6 of this report. The resulting GLS-derived regression equations follow the general structure shown below:

$Q_P=a{A}^b{B}^c$

(1)

where

Q_p

is the dependent variable, the P-percent AEP, in cubic feet per second (ft³/s);

A and B

are independent (explanatory) variables; and

a,b, and c

are regression coefficients.

If the dependent variable, Q_p, and the independent variables, A and B, are logarithmically transformed, the model takes the following form:

$\log Q_P=\log \left(a\right)+b(\log \left(A\right))+c(\log \left(B\right))$

(2)

GLS-derived regression equations for Coastal Plain and Piedmont subregions:Coastal Plain

$Q_{50\text{\%}}=\left({10}^{1.526}\right)\left(DRNARE{A}^{0.723}\right)\left(BSLDEM3M^{0.981}\right)\left({\left(SOILA+1\right)}^{-0.257}\right)$

(3)

$Q_{20\text{\%}}=\left({10}^{1.932}\right)\left(DRNARE{A}^{0.708}\right)\left(BSLDEM3M^{0.774}\right)\left({\left(SOILA+1\right)}^{-0.299}\right)$

(4)

$Q_{10\text{\%}}=\left({10}^{2.138}\right)\left(DRNARE{A}^{0.702}\right)\left(BSLDEM3M^{0.679}\right)\left({\left(SOILA+1\right)}^{-0.318}\right)$

(5)

$Q_{4\text{\%}}=\left({10}^{2.355}\right)\left(DRNARE{A}^{0.697}\right)\left(BSLDEM3M^{0.586}\right)\left({\left(SOILA+1\right)}^{-0.337}\right)$

(6)

$Q_{2\text{\%}}=\left({10}^{2.492}\right)\left(DRNARE{A}^{0.695}\right)\left(BSLDEM3M^{0.531}\right)\left({\left(SOILA+1\right)}^{-0.347}\right)$

(7)

$Q_{1\text{\%}}=\left({10}^{2.616}\right)\left(DRNARE{A}^{0.693}\right)\left(BSLDEM3M^{0.483}\right)\left({\left(SOILA+1\right)}^{-0.356}\right)$

(8)

$Q_{0.5\text{\%}}=\left({10}^{2.729}\right)\left(DRNARE{A}^{0.691}\right)\left(BSLDEM3M^{0.440}\right)\left({\left(SOILA+1\right)}^{-0.364}\right)$

(9)

$Q_{0.2\text{\%}}=\left({10}^{2.866}\right)\left(DRNARE{A}^{0.689}\right)\left(BSLDEM3M^{0.389}\right)\left({\left(SOILA+1\right)}^{-0.374}\right)$

(10)

Piedmont

$Q_{50\text{\%}}=\left({10}^{2.777}\right)\left(DRNARE{A}^{0.750}\right)\left(BSLDEM3M^{-0.679}\right)\left({\left(IMPERV+1\right)}^{0.216}\right)\left({\left(STORNHD+1\right)}^{-0.665}\right)$

(11)

$Q_{20\text{\%}}=\left({10}^{3.014}\right)\left(DRNARE{A}^{0.724}\right)\left(BSLDEM3M^{-0.609}\right)\left({\left(IMPERV+1\right)}^{0.166}\right)\left({\left(STORNHD+1\right)}^{-0.690}\right)$

(12)

$Q_{10\text{\%}}=\left({10}^{3.140}\right)\left(DRNARE{A}^{0.711}\right)\left(BSLDEM3M^{-0.570}\right)\left({\left(IMPERV+1\right)}^{0.139}\right)\left({\left(STORNHD+1\right)}^{-0.702}\right)$

(13)

$Q_{4\text{\%}}=\left({10}^{3.278}\right)\left(DRNARE{A}^{0.698}\right)\left(BSLDEM3M^{-0.529}\right)\left({\left(IMPERV+1\right)}^{0.110}\right)\left({\left(STORNHD+1\right)}^{-0.715}\right)$

(14)

$Q_{2\text{\%}}=\left({10}^{2.658}\right)\left(DRNARE{A}^{0.657}\right)\left({\left(SOILD+1\right)}^{0.188}\right)\left({\left(IMPERV+1\right)}^{0.195}\right)\left({\left(STORNHD+1\right)}^{-0.814}\right)$

(15)

$Q_{1\text{\%}}=\left({10}^{2.746}\right)\left(DRNARE{A}^{0.650}\right)\left({\left(SOILD+1\right)}^{0.200}\right)\left({\left(IMPERV+1\right)}^{0.180}\right)\left({\left(STORNHD+1\right)}^{-0.824}\right)$

(16)

$Q_{0.5\text{\%}}=\left({10}^{2.826}\right)\left(DRNARE{A}^{0.645}\right)\left({\left(SOILD+1\right)}^{0.212}\right)\left({\left(IMPERV+1\right)}^{0.166}\right)\left({\left(STORNHD+1\right)}^{-0.834}\right)$

(17)

$Q_{0.2\text{\%}}=\left({10}^{2.922}\right)\left(DRNARE{A}^{0.638}\right)\left({\left(SOILD+1\right)}^{0.227}\right)\left({\left(IMPERV+1\right)}^{0.149}\right)\left({\left(STORNHD+1\right)}^{-0.848}\right)$

(18)

Table 3.    

Number of streamflow-gage sites included in the peak-flow regression analyses by subregion and state.

[POR, period of record; ≥, greater than or equal to; DE, Delaware; MD, Maryland; NJ, New Jersey; PA, Pennsylvania]

Table 3.    Number of streamflow-gage sites included in the peak-flow regression analyses by subregion and state.

Region		State
		DE	MD	NJ	PA	Total
Coastal Plain	Total stations	22	16	16	0	54
	POR≥25	9	7	10	0	26
Piedmont	Total stations	7	8	0	25	40
	POR≥25	3	4	0	17	24
All	Total stations	29	24	16	25	94
	POR≥25	12	11	10	17	50

Table 4.    

Summary of drainage area, number of streamflow-gage sites, and average years of record used in the peak-flow regression analyses for Delaware.

[--, not applicable]

Table 4.    Summary of drainage area, number of streamflow-gage sites, and average years of record used in the peak-flow regression analyses for Delaware.

Drainage area (square miles)	Number of streamflow-gage sites			Average years of observed data record
	Coastal Plain	Piedmont	Total	Coastal Plain	Piedmont	Total
0–1	3	2	5	13.0	32.7	22.8
1–2	3	3	6	23.6	12.0	17.8
2–5	10	3	13	26.7	30.0	27.5
5–10	16	8	24	21.8	28.0	23.8
10–20	9	5	14	38.1	34.8	36.9
20–50	8	11	19	49.4	30.5	38.4
50–100	4	6	10	47.5	72.6	63.1
100–200	1	1	2	84.0	18.0	51.0
200–500	--	1	1	--	71.0	71.0
Total	54	40	94	31.8	35.2	33.2

Table 5.    

Ranges of basin characteristics used to develop the peak-flow regression equations.

[--, data not used in equations; NHD, National Hydrography Dataset]

Table 5.    Ranges of basin characteristics used to develop the peak-flow regression equations.

Basin characteristic	Piedmont		Coastal Plain
	Minimum	Maximum	Minimum	Maximum
Drainage area (square miles)	0.4	318.5	0.5	112.3
Mean basin slope (percent)	4.7	17.4	1.4	6.5
Hydrologic soil type A (percent)	--	--	0	82.4
Hydrologic soil type D (percent)	0	38.8	--	--
Storage (NHD) (percent)	0	6.0	--	--
Impervious area (percent)	0.7	40.2	--	--

Table 6.    

Values needed to determine 90-percent prediction intervals for the best peak-flow regression equations by subregion in Delaware.

[t, critical value from Student’s t-distribution for the 90-percent probability; MEV, regression model-error variance; U, covariance matrix of model coefficients; Intercept, y-axis intercept of regression equation; --, not applicable]

Table 6.    Values needed to determine 90-percent prediction intervals for the best peak-flow regression equations by subregion in Delaware.

Annual-exceedance probability	Coastal Plain							Piedmont
	t	MEV	U	U	U	U	U	t	MEV	U	U	U	U	U	U
50	1.676	0.046572	--	DRNAREA	BSLDEM3M	SOILA	Intercept	1.690	0.011247	--	DRNAREA	BSLDEM3M	IMPERV	STORNHD	Intercept
50	--	--	DRNAREA	4.14E-03	1.32E-03	–0.00062	–0.00422	--	--	DRNAREA	1.25E-03	–0.0011	–6.22E-05	–0.00125	–0.00042
50	--	--	BSLDEM3M	1.32E-03	4.87E-02	0.002959	–0.02837	--	--	BSLDEM3M	–1.10E-03	0.033579	2.94E-03	–0.00195	–0.03192
50	--	--	SOILA	–6.20E-04	2.96E-03	0.004968	–0.00616	--	--	IMPERV	–6.22E-05	0.002942	2.30E-03	0.000686	–0.00484
50	--	--	Intercept	–4.22E-03	–2.84E-02	–0.00616	0.026367	--	--	STORNHD	–1.25E-03	–0.00195	6.86E-04	0.017082	–0.0003
50	--	--	--	--	--	--	--	--	--	Intercept	–4.22E-04	–0.03192	–4.84E-03	–0.0003	0.035416
20	1.676	0.041077	--	DRNAREA	BSLDEM3M	SOILA	Intercept	1.690	0.013727	--	DRNAREA	BSLDEM3M	IMPERV	STORNHD	Intercept
20	--	--	DRNAREA	0.003863	0.00149	–0.00057	–0.00417	--	--	DRNAREA	1.60E-03	–0.00148	–6.65E-05	–0.00152	–0.00049
20	--	--	BSLDEM3M	0.00149	0.046685	0.002551	–0.02732	--	--	BSLDEM3M	–1.48E-03	0.043803	3.76E-03	–0.00241	–0.04158
20	--	--	SOILA	–0.00057	0.002551	0.004486	–0.0055	--	--	IMPERV	–6.65E-05	0.003759	2.97E-03	0.000943	–0.00627
20	--	--	Intercept	–0.00417	–0.02732	–0.0055	0.025352	--	--	STORNHD	–1.52E-03	–0.00241	9.43E-04	0.021379	–0.00059
20	--	--	--	--	--	--	--	--	--	Intercept	–4.93E-04	–0.04158	–6.27E-03	–0.00059	0.046167
10	1.676	0.03893	--	DRNAREA	BSLDEM3M	SOILA	Intercept	1.690	0.015416	--	DRNAREA	BSLDEM3M	IMPERV	STORNHD	Intercept
10	--	--	DRNAREA	0.004006	0.00169	–0.00059	–0.00445	--	--	DRNAREA	1.89E-03	–0.00185	–7.84E-05	–0.00172	–0.00052
10	--	--	BSLDEM3M	0.00169	0.049527	0.00249	–0.02895	--	--	BSLDEM3M	–1.85E-03	0.052614	4.43E-03	–0.00273	–0.04986
10	--	--	SOILA	–0.00059	0.00249	0.00458	–0.00556	--	--	IMPERV	–7.84E-05	0.004434	3.55E-03	0.001172	–0.00746
10	--	--	Intercept	–0.00445	–0.02895	–0.00556	0.02677	--	--	STORNHD	–1.72E-03	–0.00273	1.17E-03	0.024756	–0.00087
10	--	--	--	--	--	--	--	--	--	Intercept	–5.18E-04	–0.04986	–7.46E-03	–0.00087	0.055336
4	1.676	0.040292	--	DRNAREA	BSLDEM3M	SOILA	Intercept	1.690	0.018115	--	DRNAREA	BSLDEM3M	IMPERV	STORNHD	Intercept
4	--	--	DRNAREA	0.004493	0.002024	–0.00066	–0.0051	--	--	DRNAREA	2.33E-03	–0.00238	–9.99E-05	–0.00206	–0.00056
4	--	--	BSLDEM3M	0.002024	0.056821	0.002613	–0.03308	--	--	BSLDEM3M	–2.38E-03	0.066222	5.49E-03	–0.00324	–0.06268
4	--	--	SOILA	v0.00066	0.002613	0.005087	–0.00609	--	--	IMPERV	–9.99E-05	0.00549	4.43E-03	0.001522	–0.0093
4	--	--	Intercept	–0.0051	–0.03308	–0.00609	0.030451	--	--	STORNHD	–2.06E-03	–0.00324	1.52E-03	0.030059	–0.00127
4	--	--	--	--	--	--	--	--	--	Intercept	–5.59E-04	–0.06268	–9.30E-03	–0.00127	0.069512
2	1.676	0.043251	--	DRNAREA	BSLDEM3M	SOILA	Intercept	1.690	0.018932	--	DRNAREA	SOILD	IMPERV	STORNHD	Intercept
2	--	--	DRNAREA	0.00496	0.002305	–0.00073	–0.00569	--	--	DRNAREA	0.002537	–0.00091	–0.00021	–0.00198	–0.00202
2	--	--	BSLDEM3M	0.002305	0.063528	0.002762	–0.03686	--	--	SOILD	–0.00091	0.007793	0.002424	–0.0032	–0.0088
2	--	--	SOILA	–0.00073	0.002762	0.005594	–0.00664	--	--	IMPERV	–0.00021	0.002424	0.005132	0.000975	–0.00723
2	--	--	Intercept	–0.00569	–0.03686	–0.00664	0.033835	--	--	STORNHD	–0.00198	–0.0032	0.000975	0.03356	–0.00112
2	--	--	--	--	--	--	--	--	--	Intercept	–0.00202	-0.0088	–0.00723	–0.00112	0.021196
1	1.676	0.047262	--	DRNAREA	BSLDEM3M	SOILA	Intercept	1.690	0.020859	--	DRNAREA	SOILD	IMPERV	STORNHD	Intercept
1	--	--	DRNAREA	0.005601	0.002643	–0.00083	–0.00646	--	--	DRNAREA	0.002866	–0.00104	–0.00024	–0.00219	–0.00229
1	--	--	BSLDEM3M	0.002643	0.072233	0.003033	–0.0418	--	--	SOILD	–0.00104	0.008898	0.00279	–0.00364	–0.01004
1	--	--	SOILA	–0.00083	0.003033	0.006313	–0.00744	--	--	IMPERV	–0.00024	0.00279	0.005867	0.001125	–0.00828
1	--	--	Intercept	–0.00646	–0.0418	–0.00744	0.038294	--	--	STORNHD	–0.00219	–0.00364	0.001125	0.037692	–0.00129
1	--	--	--	--	--	--	--	--	--	Intercept	–0.00229	–0.01004	–0.00828	–0.00129	0.024191
0.5	1.676	0.052952	--	DRNAREA	BSLDEM3M	SOILA	Intercept	1.690	0.022874	--	DRNAREA	SOILD	IMPERV	STORNHD	Intercept
0.5	--	--	DRNAREA	0.0063	0.003006	–0.00093	–0.00729	--	--	DRNAREA	0.003207	–0.00117	–0.00028	–0.00242	–0.00256
0.5	--	--	BSLDEM3M	0.003006	0.081653	0.003335	–0.04714	--	--	SOILD	–0.00117	0.010052	0.003175	–0.0041	–0.01134
0.5	--	--	SOILA	–0.00093	0.003335	0.007102	–0.00834	--	--	IMPERV	–0.00028	0.003175	0.006632	0.001279	–0.00938
0.5	--	--	Intercept	–0.00729	–0.04714	–0.00834	0.043126	--	--	STORNHD	–0.00242	–0.0041	0.001279	0.041997	–0.00146
0.5	--	--	--	--	--	--	--	--	--	Intercept	–0.00256	–0.01134	–0.00938	–0.00146	0.027311
0.2	1.676	0.061914	--	DRNAREA	BSLDEM3M	SOILA	Intercept	1.690	0.026068	--	DRNAREA	SOILD	IMPERV	STORNHD	Intercept
0.2	--	--	DRNAREA	0.007354	0.003534	–0.00108	–0.00852	--	--	DRNAREA	0.003714	–0.00136	–0.00032	–0.00278	–0.00296
0.2	--	--	BSLDEM3M	0.003534	0.095613	0.003817	–0.05508	--	--	SOILD	–0.00136	0.011762	0.003745	–0.00478	–0.01328
0.2	--	--	SOILA	–0.00108	0.003817	0.0083	–0.0097	--	--	IMPERV	–0.00032	0.003745	0.007764	0.001503	–0.01101
0.2	--	--	Intercept	–0.00852	–0.05508	–0.0097	0.050321	--	--	STORNHD	–0.00278	–0.00478	0.001503	0.048517	–0.0017
0.2	--	--	--	--	--	--	--	--	--	Intercept	–0.00296	–0.01328	–0.01101	–0.0017	0.03191

Table 7.    

Mean and median percent differences between peak-flow magnitudes at defined frequencies computed from the annual peak-flow data for streamflow-gage sites included in both the current study and previous studies.

Table 7.    Mean and median percent differences between peak-flow magnitudes at defined frequencies computed from the annual peak-flow data for streamflow-gage sites included in both the current study and previous studies.

Annual-exceedance probability	All sites		Coastal Plain		Piedmont
	Mean percent difference	Median percent difference	Mean percent difference	Median percent difference	Mean percent difference	Median percent difference
All sites
50	–0.9	0.0	–2.4	0.1	1.1	–0.2
20	5.3	0.5	6.9	0.9	3.2	0.2
10	9.2	2.0	12.7	2.4	4.6	1.2
4	13.8	3.7	19.4	5.3	6.4	2.7
2	17.2	4.6	24.4	6.9	7.7	3.8
1	20.4	6.0	29.1	8.2	9.1	4.6
0.5	23.7	7.9	33.7	10.8	10.4	6.0
0.2	27.9	9.3	39.9	13.2	12.2	7.6
Sites with no new data since last study
50	1.8	–0.5	3.3	0.0	0.1	–0.9
20	8.5	0.1	12.8	0.5	3.2	–0.1
10	12.5	1.0	18.5	2.1	5.2	0.9
4	17.1	2.5	24.6	2.4	7.8	2.6
2	20.3	3.4	28.9	3.2	9.6	3.5
1	23.3	4.6	32.9	4.5	11.5	4.7
0.5	26.2	6.3	36.6	8.3	13.3	6.2
0.2	29.9	8.3	41.4	8.7	15.8	7.9
Sites with new data since last study
50	–3.7	3.1	–7.9	0.7	0.2	3.4
20	2.0	1.9	1.1	1.5	0.3	3.4
10	5.8	4.2	7.1	4.3	0.3	3.3
4	10.5	7.1	14.4	7.6	0.4	4.2
2	14.0	7.5	20.0	10.6	0.4	4.3
1	17.5	8.9	25.5	13.2	0.5	4.5
0.5	21.1	11.1	30.9	16.1	0.5	4.7
0.2	25.9	13.2	38.4	19.6	0.5	5.0

Comparison of Results with Previous Study

Median percent differences in peak-flow magnitudes computed from the annual peak-flow data from streamflow-gage sites included in this study show that flow estimates were generally within 10 percent of those from Ries and Dillow (2006), with slightly higher estimates for the sites and years included in this report (table 7). Sites in the Coastal Plain appeared to have higher median percent differences, particularly sites with additional recorded data since the prior study.

In Ries and Dillow (2006), variables to predict peak-flow magnitudes for the Piedmont region included drainage area, the percent area forested, the percent area with well drained soils (soil A), the percent area with impervious cover, and the percent area of surface water storage. For the Coastal Plain, variables used included drainage area, mean basin slope, and percent area with well drained soils (soil A). Similarly, in this report, for the Piedmont drainage area, mean basin slope, percent area of surface water storage, and percent area with impervious cover were chosen, but percent area poorly drained soils (soil D) was selected instead of well drained soils (soil A) and percent area forested. Variables chosen for the Coastal Plain were identical to the prior report. The model errors for the equations generated in this report (table 8) were generally within 5 percent of those from the prior report, with slightly higher model error in the Piedmont subregion and slightly lower error for the Coastal Plain subregion.

Accuracy and Limitations of Peak-Flow Regression Equations

The accuracy of a regression equation depends on the model error and the sampling error associated with the data used in its development. Model error measures the ability of a set of explanatory variables to estimate the values of peak-flow characteristics calculated from the recorded streamflow-gage data used to develop the equation (Ries and Dillow, 2006). Sampling error measures the ability of a finite number of sites with a finite number of recorded annual peak flows to describe the true peak-flow characteristics of a site. Model error depends on the number and predictive power of the explanatory variables in a regression equation. Sampling error depends on the number of sites and the amount of peak-flow data at each site used in the analysis. Sampling error decreases as the number of sites and amount of data increase (Ries and Dillow, 2006). By definition, approximately two-thirds of the flow estimates derived using a given equation will contain errors within the standard error (estimate or prediction) associated with that equation.

Traditional measures of the accuracy of peak-flow regression equations are the standard error of estimate and the standard error of prediction. The standard error of estimate is derived from the model error, and is a measure of how well the estimated peak discharges generated using a regression equation agree with the peak-flow statistics generated from the at-site annual peak-flow data used to create the equation (Ries and Dillow, 2006). The standard error of prediction is derived from the sum of the model error and sampling error and is a measure of how accurately the estimated peak discharges generated using an equation will be able to predict the true value of peak discharge at the selected AEP (Ries and Dillow, 2006). Average standard errors of estimate and prediction associated with equations 3 through 18 are shown in table 8.

For the Coastal Plain region, Thomas and Sanchez-Claros (2019) similarly developed regional regression equations for 50- to 0.2-percent annual-exceedance probabilities using Bulletin 17B guidance and data through 2017 (Interagency Advisory Committee on Water Data, 1982). The study area overlapped well with sites used in the present study, including only sites in the Coastal Plain. Their equations included drainage area, watershed slope, and percent area in soil class A, which were the same variables used by Coastal Plain equations in the present study, and had similar values of standard error as equations used in this report. For the Piedmont region, Thomas and Moglen (2016) similarly developed regional regression equations for 50- to 0.2-percent annual-exceedance probabilities using Bulletin 17B guidance and data through 2012 (Interagency Advisory Committee on Water Data, 1982). The study area only partly overlapped with the sites in the present study including sites in the Appalachian Plateau, Allegheny Ridge, Blue Ridge, and Great Valley regions in addition to the Piedmont. Drainage area, limestone as a percent of watershed area, impervious area, and forest as percent of watershed areas were included in the regional regression equations—and though the present report used drainage area, watershed slope, percent watershed in soil class D, impervious area and storage along the NHD network—standard error for both sets of regional regression equations were similar.

A measure of uncertainty that can provide information about the range of error associated with an estimate is the prediction interval. The prediction interval states the minimum and maximum values within which there is a defined probability that the true value of an estimated variable exists. For instance, defining the 90-percent prediction interval for the 1-percent-chance AEP at an ungaged site would allow the user to know the minimum and maximum flow magnitudes that have a 90-percent chance of bounding the range in which the true flow value exists. A full explanation, including example computations, of the prediction interval can be found in Ries and Dillow (2006; p. 24–28).

The regression equations can be used to estimate peak-flow magnitudes for ungaged sites with natural flow conditions in Delaware. The equations should not be applied to streams with substantial flood-retention storage upstream from sites of interest (Ries and Dillow, 2006).

Additionally, in this report, static impervious area from the National Land Cover Database [NLCD] (2016) dataset was used to represent the role of urban effects in the Piedmont regression equations (Homer and others, 2012). Although the NLCD dataset has temporal coverage back to the early 2000s, impervious percent estimates at high spatial resolution are more uncertain before that time. In an effort to assess the sensitivity of peak flow estimates to changes in impervious area through time, we artificially increased the impervious percent of each watershed in each Piedmont regression equation to examine the changes in peak-flow estimates. Doubling the impervious area of a watershed on average resulted in a 7–14 percent increase in the flow estimate depending on the AEP used (table 9). Future studies may utilize alternate datasets with proxy information on impervious area going back many decades (for example, Uhl and others, 2021) to provide an estimate of urban effects on peak flows through time.

The accuracy of the regression equations is known only within the range of basin characteristics that were used to develop the equations. The equations can be applied to ungaged sites with basin characteristics that are not within the range of applicability, but the accuracy of the resulting estimates will be unknown (Ries and Dillow, 2006).

Procedures for Estimating Peak-Flow Magnitudes

At ungaged sites on ungaged streams with no current or historical recorded at-site data, the best estimates of peak-flow characteristics can be obtained from regression equations 3 through 18 presented in this report. At gaged sites, and ungaged sites near gaged sites on the same stream, the best estimates of peak-flow characteristics can be obtained by using the results of those equations in conjunction with recorded data from a streamflow-gage site as described below.

The best magnitude estimates for peak flows are often obtained through a weighted combination of estimates produced from more than one method (Ries and Dillow, 2006). Tasker (1975) demonstrated that if two independent estimates of a streamflow statistic are available, a properly weighted average of the independent estimates will provide an estimate that is more accurate than either of the two independent estimates. Improved magnitude estimates for peak flows can be determined for Delaware streamflow-gage sites by weighting estimates from the at-site annual peak-flow data from each site with estimates obtained from the regression equations provided in this report. Improved estimates can be determined for ungaged sites in Delaware by weighting the estimates obtained from the regression equations with estimates determined based on the flow per unit area of an upstream or downstream streamflow-gage site, when available (Ries and Dillow, 2006).

In this section, step-by-step examples are provided that demonstrate how to compute magnitude estimates for peak flows. The examples are based on the methods described previously in the section “Estimating Peak-Flow Magnitudes for Ungaged Sites.” Basin characteristics that are needed for use in the estimation equations can be determined as described in table 2.

Determination of Streamflow Analysis Regions

Similar to the peak-flow analysis, prior reports (Carpenter and Hayes, 1996; Doheny and Banks, 2010) demonstrated that separating low-flow sites into Coastal Plain and Piedmont physiographic regions resulted in improved regression equation performance. In this study, we implement the same separation and note improved regression performance when data are separated by physiographic region as described in the regression equation development section below.

Flow Statistics for Streamflow-Gage Sites

Sixty-eight streamflow-gage sites were identified within the State of Delaware, or within a 25-mi buffer of the Delaware border (located within Maryland, New Jersey, or Pennsylvania) that each had 10 or more complete climatic years (April 1 to March 30) of recorded daily flow data (fig. 5, table 10). Climatic years were used in line with common practice in prior low-flow studies (Carpenter and Hayes, 1996; Dudley and others, 2020) as this period better brackets the occurrence and duration of low-flow periods than the calendar year or water year. No partial-record streamflow-gage sites were used owing to the uncertainty associated with implementing correlations in flow between sites, and no sites containing minimum flows of zero were included in the analysis. The 7-, 14-, and 30-consecutive-day low-flow discharges for annual non-exceedance probabilities (ANPs) of 50, 10, and 5 percent (corresponding to the 2-, 10-, and 20-year recurrence intervals, respectively) were calculated using the USGS SW Toolbox (Kiang and others, 2018). Forty-five of these streamflow-gage sites were subsequently used to develop low-flow regression equations following analyses of redundancy and regulation as described in the sections below. Note that although several sites included in the analyses exhibit infrequent and intermittent periods of recorded stage data that are affected by tidal influences, the published daily discharge data associated with those sites (and used in the analyses) describe computed flows (and flow statistics) that are not affected by the tidally influenced stage data. Flow statistics are provided for sites included in the regression analyses as well as those not suitable for inclusion as a result of regulation or redundancy. At-site, regression-based and weighted flow statistics are provided for sites included in regression equation development, whereas only at-site statistics are provided for sites not included in regression equation development. Weighted flow statistics were calculated as in the section “Procedures for Weighting with Regional Regression Equations.” At-site, regression-based and weighted estimates of low-flow statistics, along with basin characteristics, can be found in “Low flow estimates” and “Low basin stats” in Hammond (2021a).

Development of Regression Equations