Purpose and Scope

This study updates the previously developed statistics and equations for estimating the magnitude and frequency of flood streamflows on rural, unregulated streams in and near Virginia and West Virginia. This study supersedes previous studies by Wiley and Atkins (2010) and Austin and others (2011) and is the first study to combine Virginia and West Virginia as a study unit. This report presents revised regional skew coefficients developed using Bayesian generalized least squares (B–GLS) regression of at-site skews developed using the expected moments algorithm (EMA) as described in England and others (2018). This report and an accompanying data release (Duda and others, 2025) document the information used to estimate flood-frequency streamflows and includes a list of regional floods in Virginia and West Virginia; a discussion of climatic conditions and drainage basin characteristics affecting flooding; a presentation of results from previous studies; and an inventory of data sources containing peak streamflows, basin characteristics, and skew coefficients. At-site statistics are presented for gaged streams, both regulated and unregulated. At-site statistics are also presented as equations for estimating the streamflows of the 0.5, 0.2, 0.1, 0.04, 0.02, 0.01, 0.005, and 0.002 AEPs (2-, 5-, 10-, 25-, 50-, 100-, 200-, and 500-year floods) on rural, unregulated streams in and near Virginia and West Virginia. The statistical methods described, including an accounting of error, provide users with the associated uncertainty of flood streamflow estimates determined from applying flood-frequency streamflow equations to Virginia and West Virginia streams. Changes since the previous studies, procedures for estimating flood-frequency streamflows, and the Virginia and West Virginia StreamStats web applications are discussed in this report.

Peak Streamflow Records

The annual peak streamflow series for streamgages are documented by U.S. Geological Survey (USGS) in the peak streamflow file in the National Water Information System (NWIS) database (U.S. Geological Survey, 2023a). Flood-frequency analysis incorporates information about the annual peak streamflow series, a time series of the instantaneous maximum streamflow recorded during a water year, for streamgages of interest in an area. Streamflow information includes both systematic and historical records. Systematic records are collected at an active streamgage. Historical records are obtained from reliable documentation of flood elevations and timing at or near streamgage locations outside of their systematic period of record.

Some streamgages were operated in Virginia and West Virginia beginning in 1876 (Grover and others, 1925; Follansbee, 1994). Statewide streamgaging networks intended for long-term data collection and supported by cooperative State and Federal funding were established in Virginia in 1925 and in West Virginia in 1930 (Follansbee, 1953). The earliest streamgages relied on local observers to make daily or twice-daily water-level measurements. During the 1920s and 1930s, continuous-record streamgages with stage recorders replaced observers or were established at new locations (Grover and others, 1925; Follansbee, 1994). The number of streamgages, and thus recorded annual peaks available for this study, increased steadily until 1965 (figs. 1 and 2) (U.S. Geological Survey, 2023a). Crest-stage streamgaging networks on predominantly small streams were developed in Virginia during the 1950s and in West Virginia during the 1960s to augment systematic records from the continuous-record streamgaging network (Miller, 1978; Runner, 1980). Crest-stage streamgages record the maximum water level between inspections and are suitable for collected annual peak-flow records.

This crest-stage streamgaging network was largely implemented between 1965 and 1970. Among the streamgages in and near Virginia and West Virginia that were considered for regression analysis for this study, the largest number of annual peaks in any water year was 561, in 1971 (fig. 2). The number of peaks decreased from 550 in 1975 to 348 in 1985. In 1986, 372 peaks were recorded, including 7 historical peaks from the flood of November 1985. Since then, the number of available annual peaks has been stable, fluctuating between 292 in 1998 and 1999 and 348 in 1985. The number of available peaks for regression analysis each year does not match the number of active streamgages, because it includes historical peaks but excludes streamgages with fewer than 10 years of record and streamgages on regulated streams.

Generally, historical floods may be better documented on major rivers than on small streams (fig. 3). There was a total of 209 recorded historical peaks representative of 857 streamgages used in unregulated frequency analysis. Of these, only 47 peaks (22.5%) were recorded amongst the 526 smaller streamgages (61.4%) that drained less than 100 mi² (table 3 in Duda and others, 2025). Historical floods in the Tennessee River and its major tributaries were inventoried in 1961, through interviews with residents, reviews of historical documents, and reviews of floodmarks (Tennessee Valley Authority, 1961). Historical floods in the rest of the study area have been documented, although less thoroughly than in the Tennessee River Basin. Sometimes, historical floodmarks are obtained when a new streamgage is established, but often, historical floods were investigated to put a major contemporary flood in a historical context.

Historical peaks may be recorded for isolated floods that are locally memorable, but they are more often recorded for widespread floods; the peak streamflow file includes 353 historical peaks for unregulated, nonurban, and non-dam-break peaks. Of these, 158 were recorded in years with 10 or more historical peaks, and 234 were recorded in years with 5 or more historical peaks (table 1). More historical peaks (24) were recorded for 1940 than in any other year because of major floods in southern Virginia and western North Carolina caused by an unnamed hurricane.

Historical peaks were formerly considered to be any annual peaks that were collected or known from periods outside systematic data collection. However, the default settings in USGS computer program PeakFQ treat the values for peaks coded as “historical” as representing large, historically significant floods, even though the peak streamflow file contains many nonsystematic peaks that represent ordinary floods. Some historical peaks were obtained from floodmarks that were available from shortly before or after a streamgage was active rather than because this flood was extraordinarily large. Some apparently historical peaks were salvaged from erroneous systematic records obtained from problematic daily observer readings because the annual peaks had been verified by a hydrographer. Some post-flood inspections near the flooded area may find evidence of an ordinary flood; for such inspections, the peak values previously were coded as “historical” for lack of a better designation. After the development of the EMA methodology, the USGS introduced a specific designation for “opportunistic” peaks, or peaks of ordinary magnitude that were documented outside systematic data collection, for use in peak-flow analyses (Sando and McCarthy, 2018). This study was the first done in either Virginia or West Virginia since the development of the opportunistic code. During the initial flood-frequency analysis, the coding of all the historical peaks was reviewed in the context of other annual peak data for that site. Historical peaks that had been frequently exceeded were recoded as opportunistic and excluded from the flood-frequency analysis.

Major Floods

At-site analyses include all annual peaks but are sensitive to the largest flood streamflows known or recorded at a streamgage, referred to as the “peak of record.” The peaks of record often represent large regional floods. Most annual peaks in both Virginia and West Virginia were caused by continental frontal systems from February to April (table 2). Many of the largest, most widespread floods, however, have been caused by tropical cyclones in the summer and fall. The most common month for the peak of record in Virginia was August, and in West Virginia was November. Tropical cyclones accounted for 7 of the 8 storms that caused the peak of record through water year 2021 at more than 15 streamgages among the 686 streamgages from Virginia and West Virginia that were considered for unregulated frequency analysis; streamgages from adjacent states were not included in this analysis (table 3). These tropical cyclones have included storms from both the Atlantic Ocean and the Gulf of America, and storms from both water bodies have crossed the Appalachian Mountains and subsequently caused major floods. The regional floods that resulted in the largest number of peaks of record among the streamgages considered for regression analysis in this study are listed later in this section, along with other selected memorable floods. The number of peaks of record caused by a given flood is a function not only of the magnitude of the flood, but also the extent of the streamgaging network at the time of the flood and the subsequent effort to obtain historical peaks documenting the flood.

Previous Studies

Methods have been developed to estimate the frequency and magnitude of annual peak streamflows for rural, unregulated streams in and near Virginia and West Virginia beginning in 1954. Tice (1954) documented this process for the Shenandoah Valley region of Virginia, using the probability method proposed by Gumbel (1945). In the 1960s, a collection of works were published as part of a national effort to regionalize flood-frequency statistics using new methods outlined by Dalrymple (1960). These works included streams within Virginia and West Virginia that spanned the South Atlantic Slope Basins (Speer and Gamble, 1964a), the Cumberland and Tennessee River Basins (Speer and Gamble, 1964b), the Ohio River Basin, except the Cumberland and Tennessee River Basins (Speer and Gamble, 1965), and the North Atlantic Slope Basins (Tice, 1968). Statewide flood frequency and magnitude estimation methods were first developed in Virginia for the large hydrologic basins in 1969 (Miller, 1969) and for the small streams program in 1978 (Miller, 1978). In West Virginia, statewide flood frequencies were computed by Frye and Runner (1969, 1970, and 1971) using methods defined in the previous reports by Speer and Gamble (1965) and Tice (1968), and those proposed by Benson (1962a, 1962b)8.

In 1965, the Water Resources Planning Act recognized the need for the Federal Government, States, localities, and private enterprises to establish a comprehensive and coordinated approach to water resource conservation, development, and use (U.S. Water Resources Council, 1967). This act led to the establishment of the Hydrology Committee of the U.S. Water Resources Council which was assigned responsibility for developing and improving methods for flood-frequency analysis and providing a set of techniques based on the best known hydrological and statistical procedures. In 1967, the Hydrology Committee published Bulletin 15, titled “A Uniform Technique for Determining Flood Flow Frequencies.” This publication established the framework for standard methodology, analysis, and reporting that could be replicated between study areas. The ongoing improvements and refinements of these techniques led to periodic updates and revisions to these standards.

In 1976, the first update was issued in Bulletin 17, which introduced methods for dealing with outliers, historical flood information, and regional skew. In 1977, Bulletin 17 was revised and reissued as Bulletin 17A, which clarified the computation of weighted skew. In 1982, the committee—now known as the Interagency Advisory Committee on Water Data—released Bulletin 17B. This version provided improvements and new techniques, including better methods for estimating and applying regional skew, weighting skew, detecting outliers, and using conditional probability (England and others, 2018).

In Virginia, updated analyses for rural, unregulated streams were provided using the techniques identified in Bulletin 17 (Miller, 1978) and Bulletin 17B (Bisese, 1995; Austin and others, 2011). These studies incorporated the generalized skews developed for the United States in Bulletins 17 and 17B, respectively. In West Virginia, rural, unregulated streams were analyzed in accordance with Bulletin 17 (Runner, 1980), and Bulletin 17B (Wiley and others, 2000). Additionally, Atkins and others (2009) developed a generalized skew study for West Virginia using Bulletin 17B methods. The skews developed for West Virginia were used in the most recent update of rural, unregulated flood-frequency statistics and equations (Wiley and Atkins, 2010).

In 2018, the Subcommittee on Hydrology for the Advisory Committee on Water Information published the most recent methods for determining flood frequencies in Bulletin 17C. This version of standards adopts a generalized representation of flood data that incorporates interval and censored data types, introduces the application of the EMA, and improves methods for computing confidence intervals (England and others, 2018).

This study of flood magnitude and frequency in Virginia and West Virginia includes improvements over previous studies by identifying new regions, computing new regional skews for the combined study area, updating analytics based on methods introduced in Bulletin 17C, and incorporating additional years of data since the most recent studies (Wiley and Atkins, 2010; Austin and others, 2011).

Description of Study Area

The core study area was Virginia and West Virginia. To increase the number of streamgages available for analysis and decrease the uncertainty as to the representativeness of streamgages at the periphery of the core study area, streamgages at the 8-digit hydrologic unit code (HUC 8s) that bordered the core study area were considered for the study; the core study area and the adjacent HUC 8s comprise the extended study area.

The U.S. Environmental Protection Agency (EPA) delineated ecoregions, areas of similar geology, landforms, soils, vegetation, climate, land use, wildlife, and hydrology, to use as a framework for research, assessment, and management of terrestrial and aquatic resources. Landforms and geology are the primary components of physiography (Omernik, 1987). Ecoregions are often aligned with physiographic province boundaries but are mapped at a finer scale and account for more factors relevant to hydrology. A recent USGS flood-frequency regionalization study that borders Virginia used EPA Level III (L3) ecoregions as its regionalization framework (Feaster and others, 2023). The extended study area of Virginia, West Virginia, and adjacent HUC 8s (fig. 3B) includes parts of 10 EPA L3 ecoregions—Erin Drift Plain, Northern Piedmont, Piedmont, Southeastern Plains, Ridge and Valley, Central Appalachians, Western Allegheny Plateau, Blue Ridge, Southwestern Appalachians, and Middle Atlantic Coastal Plain (fig. 4) (U.S. Environmental Protection Agency, 2013). The 10 EPA L3 ecoregions are described from west to east (fig. 4). The climate of the entire study area is humid continental with hot or warm summers (Kottek and others, 2006). Precipitation and precipitation intensity are discussed by 10 L3 ecoregions.

A small part of the extended study area in eastern Kentucky is within the Southwestern Appalachians, an area of low mountains, hills, and intervening valleys (Woods and others, 2002). Moderate- to high-gradient streams are common and have cobble- or boulder-dominated substrates. Low-gradient streams also occur and have gravelly or sandy bottoms. Land cover is predominantly mixed-age forest (Woods and others, 1996). In the part of the Southwestern Appalachians within the study area, the 100-year 24-hour precipitation intensity (I100Y24H) ranges from 6.02 to 6.36 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) within the study area ranges from 53 to 54 in. (PRISM Climate Group, 2021).

Another small part of the extended study area is the Erie Drift Plain ecoregion in eastern Ohio. The ecoregion is glaciated and characterized by ground moraines and rolling terrain. Most soils are poorly drained. The climate in the area is continental and not influenced by Lake Erie (Woods and others, 1996). Within the Erie Drift Plain in the study area, the I100Y24H ranges from 4.92 to 5.10 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) within the same area ranges from 40 to 41 in. (PRISM Climate Group, 2021).

The Western Allegheny Plateau, in southwestern Pennsylvania, southeastern Ohio, western West Virginia, and northeastern Kentucky, is a low, dissected, and mostly unglaciated plateau. From west to east, the terrain gradually becomes more rugged. Maximum elevations are less than 2,000 feet (ft; North American Vertical Datum of 1988) and local relief is approximately 200 to 550 ft. It is composed of horizontally bedded sandstone, shale, siltstone, limestone, and coal of Pennsylvanian and Permian age. The potential natural vegetation is primarily Appalachian oak forest (Woods and others, 1996; Woods and others, 2002; Duda and others, 2025). Within the Western Allegheny Plateau in the study area, the I100Y24H ranges from 4.84 to 5.97 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) in the study area ranges from 40 to 58 in., with a median of 46 in. (PRISM Climate Group, 2021).

The Central Appalachians ecoregion includes parts of south-central Pennsylvania, western Maryland, central West Virginia, southwest Virginia, eastern Kentucky, and eastern Tennessee. This ecoregion is a high, dissected, and rugged plateau made up of sandstone, shale, conglomerate, and coal of Pennsylvanian and Mississippian age. The plateau is locally punctuated by a limestone valley and a few anticlinal ridges (Woods and others, 1996). Local relief varies from less than 50 ft in mountain glades to over 1,950 ft. High gradient streams are common. Maximum elevations generally increase towards the east and range from about 1,200 to 4,600 ft. Elevations can be high enough to ensure a short growing season, a great amount of rainfall, and extensive forest cover (Woods and others, 1996; Duda and others, 2025). Within the Central Appalachians in the study area, the I100Y24H ranges from 4.62 to 7.83 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) in the study area ranges from 39 to 73 in., with a median of 50 in. (PRISM Climate Group, 2021).

The Ridge and Valley ecoregion extends from New York to Alabama, and within the study area comprises large parts of central Pennsylvania, western Maryland, eastern West Virginia, western Virginia, and northeastern Tennessee. It is characterized by alternating forested ridges and agricultural valleys that are elongated and folded and faulted. Elevations range from about 500 to 4,800 ft. Local relief varies widely from approximately 50 to 1,500 ft. Valleys were formed from relatively weak rock, generally shale and limestone. Caves, springs, and sinkholes are common in the limestone valleys. In contrast with the rest of the study area, most streamgaging networks are trellised (Omernik, 1987; Woods and others, 1996; Duda and others, 2025). Within the part of the Ridge and Valley in the study area, the I100Y24H ranges from 4.16 to 8.71 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) in the study area ranges from 35 to 67 in., with a median of 46 in. (PRISM Climate Group, 2021). The minimum of both I100Y24H and mean 30-year precipitation in the entire study area are in the Ridge and Valley ecoregion. The minimum I100Y24H is in southern West Virginia in parts of the Bluestone and Greenbrier River Basins. The minimum mean 30-year precipitation is in areas of the South Branch Potomac and Shenandoah River Basins that are in a rain shadow from the mountains to their west.

The Blue Ridge Mountains ecoregion is a narrow strip of mountainous ridges that are forested and well-dissected (Omernik, 1987; Woods and others, 1996). This ecoregion runs from Pennsylvania to Georgia and, in the study area, comprises a narrow strip of southeastern Pennsylvania, central Maryland, the easternmost tip of West Virginia, and north-central Virginia. South of the Roanoke River, the Blue Ridge Mountains widen through south-central Virginia, western North Carolina, and eastern Tennessee. North of the Roanoke River, just three different rock types form the crest of the Blue Ridge Mountains. Maximum elevations are about 1,000 ft. In contrast, south of the Roanoke River, the Blue Ridge Mountains are higher and have more complex lithology. Maximum elevations are more than 5,700 ft near Mount Rogers. Local relief is high, and both the side slopes and channel gradients are steep (Woods and others, 1996; Duda and others, 2025). Within the part of the Blue Ridge in the study area, the I100Y24H ranges from 4.75 to 12.4 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) in the study area ranges from 39 to 77 in., with a median of 50 in. (PRISM Climate Group, 2021). The maximums of both the I100Y24H and mean 30-year precipitation in the study area are in the Blue Ridge ecoregion near Charlottesville, Va.

The Piedmont ecoregion is a transitional area between the mostly mountainous ecoregions of the Appalachians to the northwest and the relatively flat Coastal Plain to the southeast (Omernik, 1987; Griffith and others, 2002). This ecoregion is a complex mosaic of metamorphic and igneous rocks of Precambrian and Paleozoic age, with moderately dissected irregular plains and some hills. The soils tend to be finer textured than those in the Coastal Plain ecoregions to the south. Once largely cultivated, much of this ecoregion has reverted to pine and hardwood forests, with increasing conversion to urban and suburban land cover (Omernik, 1987; Griffith and others, 2002). Within the part of the Piedmont in the study area, the I100Y24H ranges from 7.17 to 10.8 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) in the study area ranges from 44 to 57 in., with a median of 46 in. (PRISM Climate Group, 2021).

The Northern Piedmont ecoregion, in southeastern Pennsylvania, east-central Maryland, and northern Virginia, consists of low rounded hills, irregular plains, and open valleys, and is underlain by metamorphic, igneous, and sedimentary rocks. Maximum elevations typically range from 325 ft on limestone to 1,050 ft on metamorphic rock. The climate is humid continental, characterized by cold winters and hot summers. Land cover is a mix of forest, agriculture, and some of the largest urban areas in the United States (Omernik, 1987; Woods and others, 1996; Duda and others, 2025). Within the part of the Northern Piedmont in the study area, the I100Y24H ranges from 6.70 to 10.3 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) in the study area ranges from 44 to 57 in., with a median of 46 in. (PRISM Climate Group, 2021).

The Southeastern Plains ecoregion extends from Louisiana and Mississippi to Delaware and, in the study area, comprises part of eastern North Carolina, eastern Virginia, and eastern Maryland (Omernik, 1987). This ecoregion is composed of irregular plains with a mixture of cropland, pasture, woodland, and forest. The sand, silt, and clay geology of this ecoregion contrasts with the older rocks of the Piedmont ecoregion. Altitudes and relief are greater than in the Middle Atlantic Coastal Plain ecoregion, but generally are less than in much of the Piedmont ecoregion. Streams have a relatively low gradient (as compared to the Piedmont ecoregion) with sandy bottoms. Within the study area, a natural feature, the Fall Line (fig. 4), is the western border of the Southeastern Plains ecoregion, separating it from the Piedmont and Northern Piedmont ecoregions (Griffith and others, 2002). Within the part of the Southeastern Plains in the study area, the I100Y24H ranges from 8.08 to 9.45 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) in the study area ranges from 44 to 51 in., with a median of 50 in. (PRISM Climate Group, 2021).

The Middle Atlantic Coastal Plain ecoregion extends from Georgia to New Jersey and, in the study area, includes the easternmost parts of North Carolina, areas bordering the Chesapeake Bay in Virginia and Maryland, and the Delmarva Peninsula (Omernik, 1987). This ecoregion consists of low-altitude flat plains with many swamps, marshes, and estuaries. Unconsolidated sediments underlie the low terraces, marshes, dunes, barrier islands, and beaches. Poorly drained soils are common, and the ecoregion has a mix of coarse and finer textured soils compared to the mostly coarse soils in much of the Southeastern Plains ecoregion. The Middle Atlantic Coastal Plain ecoregion is typically lower, flatter, and more poorly drained than the Southern Coastal Plain ecoregion. Streams in the Middle Atlantic Coastal Plain ecoregion have the lowest gradients and are the most sinuous in the study area (Griffith and others, 2002). Within the part of the Middle Atlantic Coastal Plain in the study area, the I100Y24H ranges from 7.98 to 11.1 in. (Bonnin and others, 2006). Mean 30-year precipitation (1991–2020 calendar years) in the study area ranges from 43 to 56 in., with a median of 50 in. (PRISM Climate Group, 2021).

Time Trends in Flood Magnitudes

The randomness of the systematic annual-peak series was statistically tested to detect a time trend using Kendall’s test for correlation, and the Kendall’s tau statistic (Kendall, 1975; Hirsch and others, 1982). Kendall’s tau and the level of significance (the probability or “p-value”) were computed with PeakFQ, version 7.5.1 (Flynn and others, 2006; Veilleux and others, 2014; U.S. Geological Survey, 2024). Kendall’s tau and the level of significance were determined for the annual peak series of 813 streamgages in the study area (table 3 in Duda and others, 2025). Streamgages were further limited to 418, each with a minimum of 30 gaged peaks, to reduce short-term effects related to streamgages that began operation in wet years and ended in dry years, or vice versa. The annual-peak series for 51 streamgages (about 12 percent of the 418 streamgages) indicated a statistically significant (p-value≤0.05) trend, with 40 of these exhibiting positive trends and 11 exhibiting negative trends. Assuming a significance level of 0.05, 21 streamgages would be expected to indicate a positive or negative trend by chance (about 5 percent of the 418 streamgages). Trends were significant at slightly more than twice that many streamgages.

Streamgages with annual peak records coded as “C” for urban were excluded from the study. However, some streamgages that drained urban areas were not coded as “urban.” A comprehensive review of urbanization coding of historical records would have required developing geospatial data to represent historical changes in urbanization, which was beyond the scope of this study. Some streamgages that drained urban areas indicated a strong positive trend, especially among smaller, more frequent peaks, that might plausibly indicate urbanization during the period of record, such as the streamgage at Chickahominy River near Providence Forge, Va. (USGS site 02042500), in the Richmond, Va. metropolitan area (fig. 5). However, many urban streamgages were established in areas that were already urbanized, where little or no trend would be expected from urbanization. Other streamgages in areas that transitioned to urbanization during the period of record indicated no trend.

Trends were compared to NLCD land-cover percentages from the National Land Cover Dataset (NLCD) (fig. 6). Of the 51 streamgages with significant trends and 30 or more years of record, 17 drained basins with 10 percent or more total developed land, as determined by the sum of the four developed land classes from the 2019 NLCD; 13 of these had positive trends. Among streamgages with 30 or more gaged peaks and 10 percent or more developed land in the basin, an ordinary least-squares (OLS) regression found a positive relationship between the value of Kendall’s tau and total development percentage (p-value less than [<] 0.001; adjusted coefficient of determination [R²]=0.162). The relationship was stronger among streamgages with 50 or more gaged peaks and 10 percent or more developed land in the basin (p-value<0.001; adjusted R²=0.550). For these 218 streamgages, 23 tau values were significant (p-value<0.05), and 22 of the significant values were positive. The greater strength of the relationship among streamgages with longer records indicates that 30 years is too short a period to reliably distinguish the effects of development from climatic cycles on trends in annual peaks at streamgages. Some of the positive trends seem likely to have been caused by something other than urbanization; about 66 percent (143 of 218) of the streamgages with 50 or more gaged peaks that drained basins with less than 10 percent developed land had a positive Kendall’s tau.

Analyses indicated that the presence of a positive trend alone was not a reliable indicator of urbanization and that few streamgages with high development and long periods of record lacked positive trends. Because the undeveloped or moderately developed streamgage basins (less than 10 percent development) included enough with significant positive trends, basins were not considered urbanized solely because of a significant positive trend. The data from these streamgages were all retained for regression analysis because of uncertainty about causation and a lack of guidance on how to account for climate change in flood frequency.

With the undeveloped and moderately developed streamgage basins excluded, the weak linear relationship between development and tau among streamgages with 50 or more gaged peaks was the strongest line of evidence. Hints of threshold effects were present among the remaining streamgage basins (fig. 7). A little of the effect was related to climate conditions throughout the periods of record (fig. 8). The threshold was not clear and distinct, and plausibly might have been any value between 20 and 40 percent development. The most recent urban flood-frequency study for Virginia found significant effects of development beginning at 10 percent, and that nearly all basins containing more than 30 percent of development indicated effects (Austin, 2014). Present analyses were broadly consistent with these patterns. Using a threshold lower than 30 percent would have excluded more sites with negative trends or ambiguous causes of trends, and there was no specific evidence that a threshold higher than 30 percent would be effective at excluding affected basins from further analysis. Streamgages draining developed areas were included in preliminary regression analysis to check if regression analysis allowed a threshold percentage of development to categorize streams as urban. These analyses are discussed later in the report, but they did not provide a clearer threshold for distinguishing the effects of development on trends in annual peaks than did the previous Virginia urban peak-flows study (Austin, 2014). All 97 streamgages with 30 percent or greater development were ultimately excluded from final regression analysis in this study, although setups and results of frequency analyses are included in the data release associated with this report (table 3 in Duda and others, 2025).

The most urbanized parts of the study area are near Washington, D.C.; Richmond, Va.; and Pittsburgh, Pennsylvania (fig. 6). Areas that have undergone urbanization were different from areas that remain undeveloped. Because of this difference, and because other choices were not available, the entire record for sites that had 30 percent or more development from the 2019 NLCD data was excluded from regression analysis. The early, pre-urbanized record presumably represented places that eventually became urban, and therefore, did not need rural regression equations to characterize present conditions. Urban regression equations are available for urbanized parts of Virginia in Austin (2014).

At-Site Frequency Analyses for Unregulated Streams

The flood-frequency streamflows at the 813 streamgages on rural, unregulated streams in and near Virginia and West Virginia were determined following the guidelines established by Bulletin 17C (England and others, 2018). The USGS computer program PeakFQ version 7.5.1 (Flynn and others, 2006) was used to perform initial frequency analyses to compute the station skew and flood-frequency streamflow estimates. Frequency analyses were developed for streamgages in and near Virginia and West Virginia (tables 3, 4, 5, and 6 in Duda and others, 2025). Many of these frequency analyses were used to determine regional skew for Kentucky, Tennessee, Virginia, and West Virginia (Wagner and others, 2025a; Wagner and others, 2025b). Determining regulation status is discussed in appendix 1.

The base-10 logarithm (log₁₀)-transformed systematic (continuous record or broken record that can be statistically treated as a continuous record) annual-peak series for each streamgage was fitted to the log-Pearson Type III probability curve. The EMA requires users to review historical peaks, perception thresholds, flow intervals, and low-outlier thresholds (England and others, 2018). Although the PeakFQ program populates some of these fields with default values, the default values are often suboptimal or inappropriate. Analytical evaluation and knowledge of streamgage history are required to appropriately fit the probability curve to the annual peak streamflow data, especially because of the importance of historical peak streamflows to frequency analysis.

Historical peaks are crucial for frequency analyses in this study because they are used to define periods with “interval data,” thereby extending a streamgage’s period of record with usable flood-related data outside the periods with systematic records. When considering historical conditions for a frequency analysis, it was assumed that the historical peak streamflows for that streamgage were not exceeded during the periods between the historical peak and the periods with systematic records. The historical values (or the values of other large floods for that streamgage) were used as the minimum “perception threshold” for times without systematic record between the historical peak and the period of systematic record (fig. 9). A “flow interval” from the perception threshold to zero was applied for those periods and incorporated in the frequency analysis. For all streamgage records used in this study, the upper perception threshold was infinity, meaning that an arbitrarily large flood would have been measured if it had occurred. For most periods at streamgages with systematic records, the lower perception threshold was zero, meaning that any flood greater than zero could have been measured, and the resulting flow interval was zero to infinity. Minimum perception thresholds were changed for cases when the lowest recordable streamflow was greater than zero, which is common for crest-stage gages. In these cases, the perceptible streamflows ranged between the minimum detectable stage and infinity. Minimum perception thresholds for crest-stage gages were not historically recorded in the NWIS database but instead were typically documented in paper station records in field offices. They were obtained for this study and used to create PeakFQ setup files and are available in the associated data release (Duda and others, 2025).

The fit of frequency curves may be affected by potentially influential low floods (PILFs). Along with the interpretation of historical peaks, identifying PILFs and their effect on the fit of the frequency curve is one of the major analytical tasks in developing flood-frequency analyses. The default method for detecting PILFs in PeakFQ is the multiple Grubbs-Beck test (England and others, 2018). Once identified, the PILFs are censored and removed from the frequency analysis. For many annual-peak series, however, the multiple Grubbs-Beck test does not exclude enough PILFs, and the resulting frequency curve fits high peaks poorly. In these cases, a single outlier threshold was identified and used to fit the frequency curve; this was done for 74 streamgages in this study (table 3, Duda and others, 2025). Common scenarios for a poor fit that can be improved by a manual PILF threshold include populations of peaks generated by mixed mechanisms, annual peak-series with one or two extraordinary floods relative to the rest of the series, crest-stage gages with poor ratings or stage record for part or all of the range of annual peaks, and crest-stage gages with a minimum recordable stage that excluded too many annual peaks. PeakFQ fails if more than half of the streamflow record is excluded, so a PILF threshold cannot be greater than the median annual peak. Annual peaks that remained within the streambanks or were below bankfull levels, however, were in many cases irrelevant to the fit of the upper tail of the frequency curve. In much of the study area, field surveys of bankfull channel characteristics suggested that between a quarter and a third of annual peaks were less than bankfull (Keaton and others, 2005; Lotspeich, 2009; Messinger, 2009). When field evidence of a bankfull condition was available, it was considered as a preliminary line of evidence when analysts fit frequency curves to poorly fitting upper tails and considered using a manual PILF threshold to improve the fit. Field evidence included (1) channel surveys made as part of studies that evaluated bankfull characteristics, including estimates of bankfull streamflow; and (2) cross sections obtained during streamflow measurements made during floods, for which a breakpoint between the channel and floodplain was taken as bankfull stage and streamflow was estimated from a rating in effect at that streamgage location and datum.

The key customization parameters for the frequency analyses are included in the PeakFQ output (table 3, Duda and others, 2025). Of the 813 streamgages with unregulated frequency analyses, 165 streamgages had 1 or more historical peaks in the analysis; 122 had 1, 34 had 2, 5 had 3, and 1 had 4. Manual PILF thresholds were specified for 74 streamgages. Frequency analyses censored PILFs for 189 streamgages, and the maximum number of PILFs at a streamgage was 60 (in Duda and others, 2025). Selected statistics representing the magnitude and variance of selected AEPs for the flood-frequency analyses for the 813 streamgages in Virginia, West Virginia, and adjacent states also are listed in table 4 in Duda and others (2025).

Regional Skew Coefficient

This section is an overview of the regional skew coefficients generated for the study area and the criteria used for the selection of streamgages to develop the regional skew. More details about the methodology used to generate the regional skew are in appendix 2.

To improve estimates of annual peak streamflows corresponding to the selected AEPs—particularly for streamgages with short annual peak-flow records (less than 30 years)—Federal guidance for flood-frequency analysis (Bulletin 17C; England and others, 2018) recommends using a weighted average of the at-site and a generalized, or regional, skew. Previous guidance (Bulletin 17B; England and others, 2018) supplied a national map of regional skew but encouraged hydrologists to develop more localized models when appropriate. Tasker and Stedinger (1986) developed a weighted least-squares (WLS) regression approach for estimating regional skew that accounted for the precision of the at-site skew, depending on the length of the record and the accuracy of a mean regional skew computed using OLS regression. Tasker and Stedinger (1989) also developed a generalized least-squares (GLS) model for hydrologic regression that accounts for the cross-correlation of streamflows among the streamgages. More recently, Reis and others (2005) introduced a Bayesian approach to GLS regression for regional analysis of skew, and Feaster and others (2009) used B–GLS regression to estimate regional skew. Veilleux (2011) presented a Bayesian weighted least squares/Bayesian generalized least squares (B–WLS/B–GLS) regression, wherein the regression parameters of regional skew models for California (Parrett and others, 2011) were determined using Bayesian weighted least squares (B–WLS) regression and the accuracy of the regression parameters and models were determined using B–GLS regression. Since 2011, B–WLS/B–GLS has been used by the USGS in regional skew studies (Veilleux and others, 2019; Veilleux and Wagner, 2019; Veilleux and Wagner, 2021; Feaster and others, 2023; Mitchell and others, 2023).

In this study, the EMA (Cohn and others, 1997) was used to estimate the at-site skew and its mean squared error of skew. The EMA allows for (1) the censoring of low annual peak streamflows by use of the multiple Grubbs-Beck test or use of manual low-outlier thresholds for screening PILFs and (2) the use of flow intervals to describe missing, censored, and historical data; therefore, the EMA complicates calculations of effective record length (and effective concurrent record length) used to describe the precision of skew estimates because the annual peak streamflows are not always represented by single values. To properly account for these complications and large cross-correlations between annual peak streamflows at pairs of streamgages, a B–WLS/B–GLS regression framework was developed to provide stable and defensible results for regional skew (Parrett and others, 2011; Veilleux, 2011). B–WLS/B–GLS uses 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, assess the model error variance and its precision, and compute various diagnostic statistics.

For this study, two unique skew regions were identified: Skew Region 1 encompasses streams in Virginia, West Virginia, and Maryland that drain to the Atlantic Ocean, including the eastern panhandle of West Virginia and most of western Maryland in appendix 2 (fig. 2.1). Skew Region 2 encompasses all the study area that drains to the Gulf of America, including Kentucky and Tennessee, most of West Virginia, and far southwestern Virginia. A total of 114 streamgages that were not redundant and had a pseudo-record length (refer to “Computing Pseudo-Record Length” in app. 2) of 50 years or greater were used to develop the final regional skew model for Skew Region 1 (refer to “VA_SkewRegion1.csv” in Wagner and others, 2025a). To explain the variability in skew, 24 basin characteristics were tested as covariates, but none provided sufficient predictive power. Therefore, a constant stable estimate of regional skew of 0.50 was selected for Skew Region 1. The average variance of prediction, 0.33 (standard error 0.574), is equivalent to the mean squared error of the regional skew and corresponds to an effective record length of 28 years (Griffis and others, 2004), an improvement over the Bulletin 17B regional skew map which has an effective record length of 17 years. A total of 387 streamgages that were not redundant and had a pseudo-record length of 30 years or greater were used to develop the final regional skew model for Skew Region 2 (refer to “VA_SkewRegion2.csv” in Wagner and others, 2025b). To explain the variability in skew, 24 basin characteristics were tested as covariates, but none provided sufficient predictive power. Therefore, a constant stable estimate of regional skew of 0.048 was selected for Skew Region 2. The average variance of prediction, 0.16 (standard error 0.400), is equivalent to the mean squared error of the regional skew and corresponds to an effective record length of 45 years, which is also an improvement over the Bulletin 17B regional skew map.

Regional Skew-Weighted Frequency Analyses for Rural, Unregulated Streams

Regional skew was used to complete frequency analyses for unregulated streams. Because of the relatively large uncertainty in the at-site skew for short to modest record lengths, the at-site skew can be combined with the regional skew to generate a weighted estimate of skew for a given streamgage basin (refer to “Weighted Skew Coefficient Estimator,” p. 25–26, and appendix 7 in England and others, 2018). All frequency analyses were recomputed using the same customization parameters described in the “At-Site Frequency Analyses for Unregulated Streams,” section, but with the key difference of applying weighted skews using either the new regional skews for Virginia, West Virginia, Kentucky, and Tennessee, or the EMA-compatible skew for the respective adjacent states (Koltun, 2019; Roland and Stuckey, 2019; Hammond and others, 202239; Feaster and others, 2023). The resulting frequency curves were compared to frequency curves produced with at-site skew, following the recommendation of England and others (2018, p. 26). A majority of the 417 long-term streamgages (those with 30 or more gaged peaks) indicated a better fit with at-site skew than with weighted skew, so the frequency curves produced using station skews were used for subsequent analyses for these streamgages. 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. Regional skew was used to weight frequency curves for the 404 streamgages with a short period of record (fewer than 30 gaged annual peaks), and these weighted frequency curves were used for subsequent analyses for these streamgages. The numbers of gaged peaks, the skew types used to produce the systematic statistics, and the systematic statistics themselves are in tables 3 and 4 in Duda and others (2025).

For streamgages in adjacent states where flood frequency and magnitude have been investigated using EMA, setup files for frequency analyses done from those states were acquired and used for this study (Koltun, 2019; Roland and Stuckey, 2019; Feaster and others, 2023). Frequency analyses for Kentucky and Tennessee are included with the frequency analyses for Virginia and West Virginia in Wagner and others (2025a and 2025b)96. Existing PeakFQ setup files were used for treatment of historical peaks, perception thresholds, flow intervals, and PILFs, and were modified only to allow the analysis of peak streamflow series prior to and including the 2021 water year. The skew options used by the original analysts were retained. The regional skews developed for the adjacent states were applied to streamgages from those states. Weighted-skew frequency analyses for streamgages on the Delmarva Peninsula used regional skew developed for Delaware (Hammond and others, 2022).

Frequency Analyses for Regulated Streams

Flood-frequency streamflows were computed for 86 streamgages on regulated streams in Virginia and West Virginia, for which 10 or more systematic peaks are available that represent present flow-regulation conditions (fig. 10; tables 5 and 6 in Duda and others, 2025). Regulated frequencies are not published for historical streamgages that were discontinued before existing major dams were built in the basins.

Streamgage record was assessed for regulation using a regulation index (RI) that was computed with a modification of the formula from Wieczorek and others (2021):

Details of the regulation assessment are discussed in appendix 1 of this report. The RI is a unitless ratio that accounts for total dam storage and drainage area in a streamgage basin. It assigns large weights to dams that have large values of both storage and drainage area relative to the streamgage drainage area, and small weights to dams that have small values of storage or drainage area (or both). Streams were considered to be regulated if the RI for a year of record exceeded 0.0040. Drainage basins were delineated and drainage areas were computed for all dams without published drainage areas in the National Inventory of Dams (NID) but with storage that would account for 1 acre-feet per square mile at any streamgage in the study (app. 1; table 7 in Duda and others, 2025). Dams without drainage areas were manually snapped to the stream grids used in the StreamStats applications for Virginia and West Virginia (U.S. Geological Survey, 2019).

Regulation conditions may change for a given streamgage as dams are built, modified, or removed in a basin. The present extent of regulation for streamflows on a stream grid from the Virginia or West Virginia StreamStats application, then delineated in StreamStats, are likely to be the most relevant condition to characterize. To accomplish this characterization, frequency analyses were done for periods of record that were homogeneous through time with the present intensity of regulation. Homogeneous periods of regulation compared to present conditions were determined primarily by reference to dam construction dates from the NID (U.S. Army Corps of Engineers, 2024). Construction dates were determined or estimated for dams that lacked construction dates in the NID for dams that accounted for 5 percent or more of the RI for streamgage record when the RI was greater than (>) 0.0040. If construction dates were available in the NID, they were used, but if the construction dates were not available in the NID, they were obtained from published documents and historical aerial photography and topographic maps (app. 1).

Once the RI had been computed and homogeneous periods of regulation had been determined, the RI was compared to the regulation coding for the peak streamflow record of all streamgages in the core study area that were downstream from dams in the NID. Codes in the peak streamflow file were corrected as necessary. Peak flow files included in the PeakFQ specification files in Duda and others (2025) were coded using the RI to identify regulated periods. The peak streamflow file in NWIS was also coded to be consistent with the peak flow files.

The computation procedure for regulated streams differed from that for unregulated streams in two principal ways: (1) at-site skew coefficients were applied to record from all streamgages without regional weighting, and (2) frequencies were computed only for homogeneous periods of regulation. Frequency curves were adjusted with manual PILF thresholds for 22 streamgages. AEP values from the frequency analyses that were computed for the regulated streamgages are available (tables 5 and 6 in Duda and others, 2025).

Streamgage Redundancy

Streamgage redundancy was determined by considering (1) the distance between basin centroids, the drainage area of each streamgage in a pair, and the drainage area ratio; (2) the record length and common years of each pair of streamgages; (3) the Spearman’s rank correlation coefficient (rho) and Pearson’s correlation coefficient (r) correlation coefficients between common, unregulated record and the percentage of overlap between pairs of basins (table 8 in Duda and others, 2025). The distances between basin centroids and the percentage of overlap between pairs of basin polygons were determined using the Near Distance and Tabulate Area geoprocessing tools, respectively, in ArcGIS 10.6 (Esri, 2017). Drainage area and record length were available from the NWIS database (U.S. Geological Survey, 2023a). Spearman’s rho and Pearson’s r were computed in R using the base and Hmisc R packages (R Core Team, 2021; Harrel and Dupont, 2023). This process of identifying redundant streamgages was adopted following discussions with the National StreamStats Coordinator (Peter M. McCarthy, U.S. Geological Survey, oral commun., 2024).

Using these metrics to identify potential redundancy, streamgage pairs were evaluated to determine which to retain for regression analysis. Pairs of streamgages with 25 percent or more overlap were screened, although when screening was complete, only members of pairs of streamgages with 50 percent or greater overlap were excluded from regression analysis. Correlation coefficients were reviewed to determine whether they provided additional usable information, but they generally confirmed patterns shown by spatial and temporal overlap. Record length and drainage area were the primary factors used to evaluate pairs of streamgages for redundancy. Longer record length was the first consideration in deciding which streamgage to keep from a pair. If the two streamgages had similar record lengths, the streamgage with a smaller basin was generally preferred to better represent smaller basins in the analysis. An exception to the basin-size criterion was that in pairs of streamgages on streams flowing into the core study area from an adjacent state, the streamgage that drained a greater percentage of the core study area was preferred. In a few cases, the secondary criteria were considered. In 2 cases when drainage basins had a high overlap percentage (>40 percent) but the streamgages had little record in common (2 years in one case and 10 years in the other), both were retained. In the Blackwater River Basin in the Virginia Coastal Plain, the streamgage at Blackwater River near Dendron, Va. (USGS site 02047500) had the most accurate stage-streamflow rating and stage record in the basin and was included in preference to other candidate streamgages.

In cases where more than two streamgages had potentially redundant relationships, they were screened as a group using the same criteria used to screen a pair of streamgages. For instance, when there were three streamgages on the same river and the middle one was potentially redundant with streamgages upstream and downstream, the advantages and disadvantages of the three were considered together. In these cases, the streamgage with more advantages was retained, even when it meant that two streamgages were excluded. Streamgages that were excluded from regression analysis for being redundant with another streamgage are identified in tables 1 and 8 in Duda and others (2025). While not included in the regression analysis, at-site flood frequencies are available for the redundant streamgages (table 4 in Duda and others, 2025). Basin characteristics were developed for streamgage basins that were not considered for this study because the geographic information system processing for extra basins required little extra effort and the characteristics are now available for other analyses. The extra streamgages included streamgages in the extended study area with either eight or nine annual peak streamflow values in the NWIS peak streamflow database and active streamgages in and near Virginia and West Virginia

Basin Characteristics

Characteristics of the basins draining to 1,050 streamgages in and near Virginia and West Virginia were developed for use as independent (or predictor) variables in regression analysis (tables 9, 10, 11, 12, 13, and 14 in Duda and others, 2025). Basin characteristics were developed if they had been significant in regression equations for flood magnitude and frequency in other parts of the eastern United States. The 121 quantitative basin characteristics included:

• characteristics of the basin polygon geometry (including drainage area, perimeter, centroid latitude and longitude), and compactness ratio;

• percentage of areas in Level III and IV ecoregions and physiographic provinces;

• climate characteristics (including 30-year mean precipitation for 1991–2020 and selected mean n-hour m-year precipitation intensity values);

• mean, maximum, and minimum elevation;

• area and percentage in each basin in each NLCD class, and sums of NLCD classes representing storage, forest, agriculture, and development;

• miles and density of hurricane tracks within basins;

• area and percentage of each basin with carbonate or unconsolidated surface geology;

• and the longest flow path and the slope of the main stream channel between points at 10 and 85 percent of its length.

Specific dataset sources and processing information are available in Duda and others (2025).

The basin characteristics for 693 streamgages (table 4; table 1 in Duda and others, 2025) were evaluated for linear correlation by using Pearson’s r and for nonlinear correlation by using Spearman’s rho (Helsel and others, 2020). The 693 streamgages were all considered candidates for regression analysis. They do not include the 111 streamgages that were identified as redundant (refer to “Streamgage Redundancy” section,) or some other streamgages that were ultimately excluded from regression analysis for various reasons discussed elsewhere in this report. The candidates for regression analysis included all 565 streamgages that were used to develop regression equations.

Basin characteristics with a logarithmic distribution were transformed using log₁₀ and evaluated for linear correlation. Both logarithmic and untransformed independent variables were evaluated in this study because some variables, such as the percentage of hydrologic regions, were not logarithmically distributed, and their transformation did not improve the linear relation with streamflow.

The set of 121 basin characteristics was screened for correlation and plausible relations with flood frequency and magnitude. An initial review of the EPA Level IV ecoregion data documented limited representation of streamgages in most of the 51 ecoregions, and they were excluded from further consideration (table 9 in Duda and others, 2025). The 70 characteristics that were retained were included in the correlation analysis; some were transformed and included as both untransformed and transformed values (table 15 in Duda and others, 2025).

Basin characteristics were evaluated for correlation using correlation and scatterplot matrices. Characteristics were considered highly correlated if the absolute value of Pearson’s r or Spearman’s rho was greater than or equal to 0.80. When one or more basin characteristics were highly correlated, the characteristic of the group with the most plausible effect on regional flood hydrology was retained for regression analysis, while the other characteristics were excluded from further analysis (table 15 in Duda and others, 2025). Three characteristics were excluded at this stage because they had implausible relations with flood magnitude. Two NLCD categories, barren land and shrub-scrub, were poorly represented in the study area and had small maximum values. Non-zero values of hurricane density had an inverse relation with flood magnitude; a low density of hurricane tracks seemed unlikely to increase rather than decrease flood magnitudes. Six other NLCD categories correlated weakly with other basin characteristics but were included in one of the four summations of NLCD categories instead of treated separately in regression analysis. The basin characteristics that were excluded from further analysis were highly correlated with one or more of the retained basin characteristics. Several groups of variables were highly cross-correlated with each other; the precipitation intensity variables were all highly correlated with each other, as were the elevation variables. Drainage area, basin perimeter, and longest flow path were highly correlated with each other. Streamgage and centroid latitude and longitude were highly correlated.

Fifteen basin characteristics were identified for evaluation in regional regression analysis (table 5). Of these, the mean 60-minute 100-year precipitation intensity and the mean 24-hour 100-year precipitation intensity were correlated at rho=0.82. Both were retained because their correlation was near the screening threshold of 0.80 and because they would be expected to affect basins of different sizes. Small basins with short times of concentration might produce peak streamflows correlated to the 60-minute 100-year precipitation intensity, whereas large basins with long times of concentration might produce peak streamflows correlated to the 24-hour 100-year precipitation intensity.

Regionalization Using Ordinary Least-Squares Regression Analysis

Regression equations and regions were defined iteratively through exploratory regression analysis. After redundant sites were excluded, 565 streamgages were retained for regression analysis. The median record length was 32 years. Several different measures of record length are available; this analysis used the number of gaged annual peak streamflows during unregulated periods of record. The median basin size was 35.8 mi² (fig. 11). Exploratory regression procedures and visualizations were performed using base R statistical packages, as well as dplyr, ggplot2, and Hmisc packages (Wickham, 2016; R Core Team, 2021; Harrel and Dupont, 2023; Wickham and others, 202398). An initial OLS regression equation was developed for the log₁₀ of the 0.01-AEP (100-year) streamflow and drainage area for all streamgages in the study area. Residuals were plotted on interactive maps, both as basin polygons and on the centroids of the basin, in ArcGIS Pro (Esri, 2024) and compared to regional boundaries from the previous cycle and 10 EPA L3 ecoregions. These ecoregions represented a better starting point than the last-cycle regions. The extended study area included parts of 10 EPA L3 ecoregions (discussed in the “Description of Study Area” section). The areas within the Erie Drift Plain and Southwestern Appalachians were small. Streamgages within the Erie Drift Plain were grouped with the Western Allegheny Plateau. Streamgages within the Southwestern Appalachians were included within the Central Appalachians for exploratory analysis but ultimately excluded from regressions.

The ecoregions were not necessarily continuous. Some residuals indicated that islands of ecoregions, exemplified by isolated pockets of the Ridge and Valley and Blue Ridge ecoregions surrounded by the Central Appalachians, Piedmont, or Northern Piedmont ecoregions, were not predictive, so regional boundaries were moved from the ecoregion border to nearby basin divides from the Watershed Boundary Dataset (WBD) (U.S. Geological Survey, 2023b). A series of regional regressions were developed for the 0.01-AEP (100-year) streamflows. They were tested in turn, and regional borders were refined to include streamgage basins with similar residuals, by adjusting the border to the HUC 8 border nearest the ecoregion border. Clusters of streamgages with similar regression residuals were grouped. Regions were reviewed on an interactive map, graphically assessed, and tested for significance as a categorical variable against the wider region. This process resulted in 12 regions (fig. 12; table 6). Many of the hydrologic region names are similar to ecoregion or physiographic province names. For clarity, the hydrologic regions are referred to by abbreviations, while the full name of ecoregions or physiographic provinces is used when discussing the ecoregion or province.

Regional boundaries were adjusted to the nearest HUC 8 borders, but there were several exceptions. The largest differences were that the entire Tygart Valley River Basin was included in Western Allegheny Plateaus (WPLT); the Cheat and Greenbrier River Basins were included in Northern Highlands (NHIGH), although large parts of their drainages are within the Ridge and Valley L3 ecoregion; parts of the Central Appalachians L3 ecoregion in the Potomac River Basin were included in Northern Valley and Ridge (NVYRD); areas of the Northern Piedmont L3 ecoregion are included in the Northern Blue Ridge (NBLU); the Delmarva Peninsula is a separate region, DLMRV; the part of the Middle Atlantic Coastal Plain ecoregion outside the Delmarva Peninsula is combined with most of the Southern Plains ecoregion in Coastal Plains (PLNS); and areas around Washington, D.C., and in northern Maryland with flood characteristics similar to the Northern Piedmont region but in the Middle Atlantic Coastal Plain were included in the Northern Piedmont (NPDMT).

The L3 ecoregion boundary between (1) the Piedmont and Northern Piedmont ecoregions and (2) the Southeast Plains ecoregion is a natural feature, the Fall Line. It was used as the boundary between the Piedmont and Costal Plains regions without adjustments. The Fall Line was an appropriate boundary because the differences in topography and stream channel slope and the sinuosity between the regions were likely to cause the differences in flood characteristics between the regions. Too many streams in this region formed in the Piedmont and drained across the Fall Line into the Coastal Plain region to make basin borders feasible as regional borders.

Although most of the regional boundaries in the Tennessee River Basin were adjusted to HUC 8 boundaries, tributaries of the Powell River in Southern Highlands (SHIGH) had flood histories and flood characteristics similar to SHIGH. The border between SHIGH and Southern Valley and Ridge (SVYRD) was adjusted by the 12-digit hydrologic unit code (HUC 12) borders so that these streams were included in SHIGH. The Bluestone River Basin drains both SVYRD and SHIGH. The regional border through the Bluestone River Basin follows the ecoregion border. Streamgages were included in the region that included the majority of their drainage area.

Three hydrologic regions that were developed from L3 ecoregions were split. The Blue Ridge was split into the NBLU and Southern Blue Ridge (SBLU) regions at the discontinuity in the Blue Ridge L3 ecoregion (figs. 4 and 12). Streamgages from the New, Tennessee, and part of the Roanoke River Basins were included in SBLU, while streamgages from the rest of the Roanoke River Basin and the James, Rappahannock, and Potomac River Basins were included in NBLU. The Ridge and Valley ecoregion was split into NVYRD, CVYRD, and SVYRD. Streamgages from the Potomac River Basin were included in NVYRD; streamgages from the James River Basin were included in the CVYRD; and streamgages from the Roanoke, New, and Tennessee River Basins were included in the SVYRD. The Central Appalachians ecoregion was split into NHIGH and SHIGH. The centerline of the New and Kanawha Rivers divides the two regions. The NHIGH includes the Cheat, Greenbrier, Gauley, and Elk River Basins, and tributaries to the New and Kanawha Rivers that enter from the right descending bank. The SHIGH includes part or all of the Bluestone, Coal, Guyandotte, Big Sandy, Kentucky, Cumberland, and Powell River Basins, and tributaries to the New and Kanawha Rivers that enter from the left descending bank.

Some streamgages in the extended study area were evaluated for the regressions and ultimately excluded because the results made the regression diagnostics worse for their prospective region, and the streamgages were on streams that drained away from the core study area. These included some streamgages in the Yadkin River Basin in North Carolina and the Cumberland River Basin in Kentucky, and all streamgages in the Youghiogheny River Basin in Maryland and Pennsylvania. Streamgages in HUC 8s adjacent to Virginia in the Cumberland and Kentucky River Basins were included in exploratory regressions. Neither basin drains any area in Virginia. Streamgage records from these HUC8s were retrieved, frequencies were computed, and streamflows were examined in preliminary regression analysis; there were no apparent graphical differences between the Cumberland and Kentucky River Basins and the rest of SHIGH. Testing these basins as a categorical variable against the rest of SHIGH found differences were not significant (p-value=0.778 and n=8 for Kentucky; p-value=0.556 and n=12 for Cumberland). The Cumberland River Basin included three streamgages in the Southern Appalachians L3 ecoregion; they were included in SHIGH for exploratory analysis but excluded from the final regressions.

For the Yadkin River Basin, three headwater streamgages that drain substantial amounts of Virginia were included in the Piedmont region. Streamgages further downstream from the North Carolina and Virginia border, with less than half of their basins in Virginia, when included, made regression diagnostics worse for their regions, and they were excluded from the final equations.

All 15 streamgages with valid unregulated frequency analyses in the Youghiogheny River Basin were excluded from the regressions for the Central Appalachian Highlands and Western Allegheny Plateaus ecoregions. The slope of an OLS regression for that specific area was significantly different from the slope of an OLS regression for the rest of the region, suggesting they would be part of a different region had the study area included more of Maryland and Pennsylvania (p-value=0.006). None of the streamgage basins had a 50 percent drainage area within West Virginia. Flood statistic estimates may be made for the part of West Virginia drained by the Youghiogheny River using the NHIGH regional equations, although with unknown errors. Users may compare the estimates for this area to estimates made using Maryland or Pennsylvania regression equations.

A broad suite of basin characteristics discussed previously was condensed to a subset of 15 basin characteristics thought to be relevant to stream flood characteristics while minimizing cross-correlation among variables. The subset of basin characteristics was tested as potential independent variables for each of the regions. None were significant as independent variables, largely because the characteristics were homogeneous within regions developed from regression residuals (fig. DR–1 of Duda and others, 2025).

The percentage of carbonate rock in the three Valley and Ridge regions was examined closely for threshold effects, which might serve as criteria for numeric valid range limits for regressions. No threshold effects were found. Streamgages with greater percentages of carbonate rock may have had higher regression residuals than streamgages draining lesser percentages of carbonate rock, but they included nearly equal numbers of positive and negative values. Streams in karst areas exchange flow with groundwater in cave and sinkhole systems; streams may either gain or lose streamflow to groundwater and affect peak streamflows. The geospatial data available to characterize carbonate rock in a consistent manner throughout the study area lacked the detail that would allow reaches that gained or lost streamflow to be analyzed separately. However, geospatial datasets that characterize karst sinkholes that were recently published or in development during this study have the potential to be valuable covariate data in future flood-frequency studies (Virginia Department of Energy, Geology and Mineral Resources Program, 2023; Shank, 2024; Ladd, 2025). These datasets, as well as geological maps that depicted sinkholes (Brezinski, 2014; Brezinski and others, 2021), were used to characterize five influential observations in Valley and Ridge regions that were ultimately excluded from regression analysis because of apparent flow loss through sinkholes (table 1 in Duda and others, 2025).

Development, as a land-cover class, was also examined closely for threshold effects, which has been discussed in the “Time Trends in Flood Magnitudes” section. Streamgages on streams that drained 30 percent or more developed land, as the sum of all development classes in the 2019 NLCD, tended to have larger regression residuals than other streamgages in the same regions, including both negative and positive residual values. Austin (2014) identified a threshold value of 30 percent development for Virginia for excluding urbanized basins from rural regressions, and this study did not find clear evidence that a different threshold was appreciably better. All streamgages at the outlets of basins that were 30 percent or more developed in 2019 were excluded from regressions (fig. 6; table 1 in Duda and others, 2025).

Streamgages with high leverage were removed from the regressions if they were unrepresentative of the regional conditions the regression equations are intended to model. Streamgages that drained more than 3,000 mi² were removed from regressions; these included main stem streamgages on the Potomac, Shenandoah, Roanoke, Monongahela, Big Sandy, Holston, James, Roanoke, Kanawha, and New Rivers. Streamgages that drained less than 0.2 mi² were all excluded because all of them had high leverage; most of them had high influence relative to their region; there were not enough of them in any single region to say whether they were representative of all small streams in the region; and other methods are available to estimate flood characteristics in streams that drain less than 0.20 mi². In general, excluding them reduced measures of variance for the regressions. The exception was the Coastal Plains region, which included a group of 4 streamgages that drained between 0.008 and 0.041 mi². As a group, they had large leverage on the regional regression equation, but when plotted, they were close to a regression line extended from streamgages that drained larger areas, indicating that this group unduly reduced measures of regression variance. They were excluded because they were unrepresentative of regional conditions; all were within a small area, the Richmond, Va. metropolitan area; drained basins more than an order of magnitude smaller than the next smallest streamgage in the region (0.99 mi²); and because they were on small streams that required different methodology than the rest of the streamgaging network (Hayes and Young, 2005).

Checks on regionalization included the examination of measured quantile compared to predicted quantile plots (plots of residuals compared to expected quantiles), plotting residuals against fitted values and the independent variable, investigating streamgages with the highest Cook’s distance values and largest residuals relative to their prospective region for errors in streamflow records and basin characteristics, checking for spatial clusters of streamgages with similar high or low residuals, and iterative review of residuals from regional regressions in interactive maps. When regions had been defined for the 0.01 AEP (100-year flood), regression equations were developed for the 0.2 AEP (5-year flood), and the same set of checks was performed for them. In every region, regression equations for the 0.2 AEP (5-year flood) had less variance than did the equation for the 0.01 AEP (100-year flood).

Development of and Accuracy of Generalized Least-Squares 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 the regional regression equations. GLS regressions were developed for the 0.5, 0.2, 0.1, 0.04, 0.02, 0.01, 0.005, and 0.002 AEPs (2-, 5-, 10-, 25-, 50-, 100-, 200-, and 500-year floods) for the 12 regions with drainage area as the only independent variable (tables 6 and 7; tables 16 and 18 in Duda and others, 2025; fig. 13). GLS regression techniques are preferred for regionalization of streamflow-frequency statistics because they improve estimates of streamflows and accuracy of the regression model by accounting for (1) unequal record lengths from streamgages and (2) cross-correlation of streamflow statistics from streamgages (Farmer and others, 2019, p. 18). Streamgages with shorter records are given less weight than streamgages with longer records. Additionally, less weight is given to streamgages where concurrent peak streamflows are correlated with nearby streamgages.

Based on results from the OLS analysis, explanatory and response variables were log-transformed prior to use in the GLS regression analysis. The weighted-multiple-linear regression (WREG) program (Farmer, 2021), 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 alpha and theta, the correlation-smoothing parameters in the weighted-multiple-linear regression program, were adjusted for each regional regression equation to improve model fit (table 6). The range of drainage areas used to develop the equations for each region is also given in table 6; the smallest valid drainage area in any region is 0.21 mi², and the largest valid drainage area in any region is 2,966 mi².

Scatterplots of observed and predicted values filtered by region were examined to assess the fit of the regression equations (fig. 14). Most regions had a good fit, although for PLNS, NPDMT, SVYRD, and possibly CVYRD, the fit was only fair and was better for larger magnitudes of observed and predicted values than for smaller magnitudes. The fit for DLMRV was poor, which was expected because of the poor performance metrics of the regression. Plots of residuals and predicted values generally indicated little unexplained structure, although the pattern of larger residuals, compared to smaller predicted streamflows, was observed for the same three regions—PLNS, NPDMT, and SVYRD—that had a worse fit for smaller than larger predicted streamflows, as did DLMRV. Measured quantile compared to predicted quantile plots were generally linear with some exceptions of for one to three points. The exceptions were NBLU, with a concave downward distribution; SVYRD and NPDMT, which indicated departures at both ends of the plot, and DLMRV, which indicated a poor fit.

Several performance metrics from the weighted-multiple-linear regression program can be used to evaluate the accuracy of the regression equations. Performance metrics for the final regression equations include the model error variance, the root mean square error, the pseudo coefficient of determination (pseudo-R²), the average variance of prediction (AVP), and the average standard error of prediction (table 8; table 18 in Duda and others, 2025). The model error variance and standard model error variance describe the portion of the total error that can be attributed to having an imperfect model, in log units and percent, respectively. The root mean square error and pseudo-R² are measures of model accuracy for streamgages used in the model development. Root mean square error describes how much the predicted peak streamflows deviate from the observed peak streamflows. 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 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 the average standard error of prediction are measures of how well the regression model performs at predicting peak streamflows for ungaged sites not used to develop the regression equations—lower values indicate greater predictive power. Pseudo-R² and AVP were the most important metrics for selecting the final regression equations. Equations for calculating these performance metrics are available in England and others (2018).

The pseudo-R² values for the final regression equations ranged from 0.481 percent (DLMRV; 0.002 AEP or 500-year flood) to 0.995 percent (SBLU; 0.5 AEP or 2-year flood) for all regions and AEPs. They ranged from 0.499 (DLMRV) to 0.968 percent (SBLU) for the 0.01-AEP (100-year) streamflows (table 8). The AVP values for the final regression equations ranged from 0.005 (SBLU; 0.5 AEP or 2-year flood) to 0.154 (NPDMT; 0.002 AEP or 500-year flood) for all regions and AEPs. They ranged from 0.026 (PDMT) to 0.115 (DLMRV) for the 0.01-AEP streamflows. The average standard error of prediction values ranged from 16.5 percent (SBLU; 0.5 AEP or 2-year flood) to 112 percent (NPDMT; 0.002 AEP or 500-year flood) for all regions and AEPs. Performance metrics were better for the statistics that characterize high-frequency floods than for those that characterize low-frequency floods; as an example, pseudo-R² values for 0.5-AEP (2-year flood) streamflows, excluding DLMRV, ranged from 0.862 (PLNS) to 0.995 (SBLU), in contrast to values for 0.002-AEP (500-year flood) streamflows, that ranged from 0.662 (PLNS) to 0.951 (SBLU). Overall, the final regression models performed poorly for DLMRV, with a pseudo-R² value of 0.4987 for the 0.01 AEP. They performed poorly to fairly well for NPDMT and PLNS, with pseudo-R² values of 0.8212 and 0.7618, respectively, and well for all other regions, which had pseudo-R² values ranging from 0.8725 to 0.9684 for the 0.01 AEP (Duda and others, 2025). Performance metrics for equations for NPDMT and PLNS show larger error terms than most of the study area and Ohio, Pennsylvania, and North Carolina (Koltun, 2019; Roland and Stuckey, 2019; Feaster and others, 2023). NPDMT and PLNS have fairly large error terms compared to many other studies, but pseudo-R² is about the same as last cycle, and standard error of prediction is somewhat to substantially larger, as discussed in the “Changes in 0.01 AEP Streamflows Since Most Recent Studies” section. These are regions with high natural variation in peak streamflows, and large error terms may indicate this natural variability. However, many streamgages in NPDMT are on streams with unstable streambanks and streambeds that make ratings unstable and difficult to develop and maintain (Messinger and Burgholzer, 2018). Because of natural rating instability, streamflow data quality throughout the region is likely to be worse than in regions where streambanks and streambeds are stable. Streambanks and streambeds in PLNS are more stable than those in NPDMT, but collecting high-quality streamflow data in the low-gradient streams typical of the region is also difficult. Station records for influential streamgages in these regions did not show specific problems in streamgage operation or the quality of indirect streamflow and other flood measurements. The basins were representative of the region according to all available basin characteristics. Arbitrarily removing the most influential streamgages from these regression analyses would have improved the regression diagnostics but would be inappropriate without a specific reason (Helsel and others, 2020).

All equations for the DLMRV are poor. The most obvious option to improve DLMRV would have been to include streamgages from Delaware and Maryland farther from the core study area in the analysis. However, Delaware has recently published updated flood-frequency regressions that include most of the streamgages that are included in the DLMRV; expanding analysis of the Delmarva Peninsula would amount to repeating that study (Hammond and others, 2022). No compelling reason was evident to repeat that study. Flood-frequency and magnitude estimates for the Delmarva Peninsula may be made with the Delaware Coastal Plain regression equations from Hammond and others (2022; app. 3), which are intended to be added to the Virginia StreamStats application (U.S. Geological Survey, 2019).

Regression-Weighted Flood-Frequency Estimates and Confidence Limits at Gaged Sites

The uncertainty in a flood-frequency estimate can be reduced by computing a weighted average of the at-site estimate and the regional regression estimate; Bulletin 17C (England and others, 2018) describes a weighting method in which the log₁₀ of the regression and at-site LP3 estimates are weighted inversely proportional to their respective variances. Q_i is the estimated peak-flow statistic at site i, and X_i is the log-transformed variable. All subsequent operations are performed on the transformed variable X_i. The weighted estimate is calculated using variances as

As described in another section of the report, “At-Site Frequency Analyses for Unregulated Streams,” PeakFQ uses the EMA to provide a direct fit of the log-Pearson Type III distribution, which includes an estimate of the variance V_site,i corresponding to each computed AEP.

For independent X_site,i and X_reg,i, the variance of the weighted estimate for site i is calculated (with all V variables in log₁₀ units) as

Confidence intervals on the weighted estimate can also be calculated. For example, upper and lower 95 percent confidence limits on the weighted quantile estimate are calculated as

and note that X_weighted,i, V_weighted,i, and CI_i must be calculated separately for each site i and for each AEP of interest.

Because of the difficulty in computing the sampling error variance, similar computations to those listed in equations 2, 3, and 4 can be used (with reduced accuracy) to compute weighted estimates and confidence limits for estimates at sites not used in model development. The only difference is that the average variance of prediction for the regression equations is substituted for V_reg. Regression-weighted estimates are provided for all streamgages where frequency analyses were computed, that were not used in valid regression analyses, estimates can be computed and for which valid regression estimates can be computed (table 19 in Duda and others, 2025).

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

A different method of combining information from gaged and ungaged streams is used for a location on an unregulated stream that is near a streamgage. If the drainage area of an ungaged site is between 50 and 150 percent of the drainage area of a gaged site on the same stream and flood-frequency characteristics have been or can be computed for the gaged site, then the following method of adjusting the estimated flood-peak streamflow of the ungaged site is available (Feaster and others, 2023; Mitchell and others, 2023):

Examples of Estimation Procedures

The techniques for estimating peak discharges described in this report can be applied to four different situations:

Examples of these procedures are described below.

2. Substitution of the basin characteristics into the equation produces:

StreamStats may be used to make estimates with the regional regression equations. If StreamStats is unavailable, the regression equations can be solved using other tools. Users who wish to use spreadsheet software may load table 16 from Duda and others (2025) in Microsoft Excel, LibreOffice Calc, or another spreadsheet, and copy the regression coefficients into formulas set up to solve the regression equations.

Limitations of Procedures for Estimating Flood-Frequency Streamflows

Procedures developed in this study are applicable only to rural, unregulated streams within the boundaries of Virginia and West Virginia. Procedures also are limited to the range of drainage areas used to develop the equations for each region, given in table 6; the smallest valid drainage area in any region is 0.21 mi² and the largest valid drainage area in any region is 2,966 mi².

Caution should be applied if flood-characteristic estimates are needed for urban areas with paved surfaces, concrete channels, or culverts that substantially alter runoff or infiltration. Streamgages used to develop regression equations drained basins that contained as much as 30 percent development from the 2019 NLCD (Dewitz and U.S. Geological Survey, 2021). However, annual peak data that were coded as affected by urbanization were excluded from the regressions and the equations published in this report are not valid for any basin that includes 30 percent or more developed area. For developed areas with between 10 and 30 percent development, users should compare results of the regression equations in this report to other available equations for estimating urban peak streamflows (Austin, 2014, for Virginia), and to investigate and consider local conditions.

The regional regression equations are not meant to be applied to streams that are regulated by dams or with large lakes and ponds that substantially retain runoff, measured by an RI greater than or equal to 0.0040. Procedures may not provide accurate estimates for extensively or intensively mined areas if excessive runoff is affected by landscape or subsurface alteration. Few streamgages were available to quantify the effects of large percentages of swamps within a basin; only two streamgages included in the PLNS region drained more than 20 percent of swamps, measured as the woody wetlands land cover class in the 2019 NLCD. Although DLMRV included 12 streamgages that drained 20 percent or more woody wetlands, the regression equations for DLMRV were poor, with pseudo-R² values ranging from 0.481 to 0.558. Regression equations from this report may not provide accurate estimates for basins that drain 20 percent or more area in the woody wetlands NLCD class.

Regression procedures may not provide accurate estimates for karst areas if excessive runoff is diverted into, outside, or deep within the basin through solution channels or other cavities in carbonate rocks. Available geospatial data depict the presence of carbonate rock, but comprehensive, quantitative information on flow paths through karst is not available. Report users who have access to detailed local information on flow paths through solution openings in carbonate rock may be able to develop localized methods of estimating flood magnitude and frequency in specific karst areas that are more accurate than regression equations in this report. Uncertainty terms for Valley and Ridge regression equations are larger than for equations for some other regions, and some of that variability may result from gains and losses to streams underlain by carbonate rock.

Equations developed for the Delmarva Peninsula for this study are poor because only one streamgage has been operated in the part of the peninsula in Virginia and only a few streamgages have been operated in the southern Delmarva Peninsula. Equations recently developed for the Coastal Plain of Delaware include the streamgages in the southern Delmarva Peninsula (Hammond and others, 2022). The equations are expected to provide better estimates of flood statistics than the equations presented in this report. These equations are included in appendix 3 of this report and are intended to be integrated in the Virginia StreamStats application.

Improved methodology may become available for applying the regression equations and weighting the at-site statistics published in this report. New methods that are deployed to the Virginia and West Virginia StreamStats web applications will have been assessed by the USGS for use in many situations but might or might not meet the needs of all users. The methods described in this report are intended to be integrated in StreamStats as the default method for estimating flood characteristics for rural, unregulated streams in Virginia and West Virginia.