Purpose and Scope

This report describes regression equations developed for estimating selected statistics for streamflows in Massachusetts from basin characteristics (hydrography, elevation, physical, land-use, soil, surficial geology, and climate). The selected streamflow statistics estimated with the regression equations are for near-natural flow conditions (minimal to no regulations). Regression equations were developed for selected streamflow statistics, including selected annual and monthly flow durations; selected monthly median flows; selected 7- and 30-day low-flow frequencies; and other statistics, such as the harmonic mean, for the mainland area of Massachusetts (fig. 1) (hereafter referred to as the “mainland” equations). Selected streamflow statistics are also provided for streamgages with regulations and streamgages with minimal to no regulations in and near Massachusetts. These statistics include the ones estimated for the regression equations and other selected annual flow durations for higher streamflows, monthly flow durations, and median of the monthly means streamflows. The streamflow statistics, basin characteristics, streamflow variability index, and regression analyses for the mainland equations are provided in a USGS data release (Bent and others, 2025).

A separate set of regression equations were developed that only estimate annual flow durations between the 50- and 99-percentiles for the Plymouth-Carver-Kingston-Duxbury aquifer system in southeastern Massachusetts and Cape Cod (fig. 1; hereafter referred to as the “southeastern” equations). Additionally, similar streamflow statistics, for which equations were not developed, are also provided for the southeastern area. The streamflow statistics, basin characteristics (aquifer, elevation, physical, land-use, soil, surficial geology, and climate), regression analyses, and other information for the southeastern equations are provided in a separate USGS data release (Carlson, 2025; Sturtevant and others, 2025).

An evaluation of the accuracies of both the mainland and southeastern equations and the limitations for their use is provided, as are considerations for further studies. Discussion about the USGS StreamStats web-based application is also provide in the report.

Previous Studies

Fennessey and Vogel (1990), Vogel and Kroll (1990), Ries (1990), Risley (1994), Ries (1994a, b, 1997, 1999), Ries and Friesz (2000), Ries and others (2000), and Archfield and others (2010) provided estimated streamflow statistics and regression equations for various flow durations such as the 7-day, 2-year low-flow frequency (7Q2) and 7Q10 in Massachusetts. These studies have included equations for low-flow frequencies and low-flow durations. Explanatory variables for the low-flow equations in these studies have included drainage area, area of stratified-drift deposits per unit of total stream length, mean basin slope, basin relief (maximum minus minimum basin elevation), average annual precipitation, open water, sand and gravel deposits, average maximum monthly temperature, X- and Y-location of basin outlet, X- and Y-location of basin centroid, and region of the State.

Wandle and Randall (1994) developed regression equations for estimating low-flow frequencies, 7Q2 and 7Q10, for high- and low-relief regions of central New England. Explanatory variables for the equations included drainage area, surficial geology, area of swamps and lakes, mean basin elevation, mean channel length, and mean annual runoff. Wandle (1983, 1987) previously developed low-flow-frequency and flow-duration equations for Massachusetts and New England, respectively. Ries (1990) developed regression equations to estimate monthly and mean annual runoff from major drainage areas in Massachusetts and Rhode Island draining to Narragansett Bay. Explanatory variables for the equations included area of till, area of stratified-drift deposits and storage (water bodies and swamps), and area of urban land. Armstrong and others (2008) provided regression equations for estimating median monthly streamflows in Massachusetts. DelSanto and others (2023) developed equations for estimating the 7Q10 for the northeastern United States (included Massachusetts streamgages) using linear regression and machine learning estimation methodologies: random forest decision trees, neural networks, and generalized additive models. Their equations included the minimum 30-day cumulative precipitation and average 30-day high temperature as well as drainage area, slope, mean elevation, wetlands area, and forest area. Bent and Archfield (2002) and Bent and Steeves (2006) provided logistic regression equations for estimating the probability of a stream flowing perennially in Massachusetts.

Regression equations for estimating selected low-flow statistics have been published in Connecticut, Rhode Island, New Hampshire, and New York studies over the last 20 years. Low-flow equations have not, currently (2025), been developed for Vermont. Ahearn and others (2006) and Ahearn (2008, 2010) and provided estimated streamflow statistics and regression equations for various flow durations, 7Q2, 7Q10, and seasonal flows based on aquatic habitat needs in Connecticut. Kliever (1996) estimated the 99-, 98-, 97-, 95-, 90-, 85-, 80-, 70-, 60-, 50-percent flow durations, 7Q10, and mean monthly streamflows for August, February, April, and May for 16 partial-record stations in northern Rhode Island. Cervione and others (1993) calculated the 99-, 98-, 95-, 90-, and 80- percent flow durations for 25 partial-record stations in southern Rhode Island. Cervione and others (1993) also presented a regression equation to estimate the 7Q10 for selected streams in Rhode Island. Bent and others (2014) provided low-flow equations for the 99- to 1-percent flow duration and the 7Q2 and 7Q10 in Rhode Island. Flynn (2003a, b) developed low-flow equations to estimate seasonal (winter, spring, summer, and fall) and annual 98-, 95-, 90-, 80-, 70-, and 60-percent flow durations and the 7Q2 and 7Q10 in New Hampshire. Randall and Freehafer (2017) developed low-flow equations for the lower Hudson River Basin, New York (area adjacent to Massachusetts and Connecticut), for the 7Q2 and 7Q10.

Description of Study Area

Low flows are greatly affected by the geography, climate, and surficial geology upstream from the measurement location. Massachusetts encompasses 8,093 mi² in the northeastern United States (fig. 1). Elevations range from sea level in coastal areas to about 3,500 feet (ft) above sea level (referenced to the North American Vertical Datum of 1988 [NAVD 88]) in the northwest. Elevations generally increase from eastern to western Massachusetts. The climate in Massachusetts is humid, with average annual precipitation ranging from about 40 to 45 inches (in.) in eastern Massachusetts to about 40 to 50 in. in western Massachusetts, where higher elevations may cause orographic effects. Average annual temperature is about 50 degrees Fahrenheit (°F) in eastern Massachusetts and about 45 °F in western Massachusetts (Bent and Waite, 2013). About half of the annual precipitation is returned to the atmosphere through evaporation and plant transpiration, with the remainder becoming groundwater recharge or stream runoff (Bent and Waite, 2013).

Surficial deposits that overlie bedrock in most of Massachusetts were deposited mainly during the last glacial period but can include areas of recent floodplain alluvium deposits along rivers and streams (Bent and Waite, 2013). In this report, these surficial deposits are classified as either till (which includes till, till with bedrock outcrops, sandy till over sand, and end-moraine deposits) or stratified deposits (which include sand and gravel, coarse sand, floodplain alluvium deposits, and fine-grained sand). Till (also known as ground moraine) is an unsorted, unstratified mixture of clay, silt, sand, gravel, cobbles, and boulders, typically deposited by glaciers on top of bedrock throughout much of the State. Till is primarily found in upland areas but can also be found at depth in river valleys. Stratified deposits include sorted and layered glaciofluvial and glaciolacustrine deposits. Glaciofluvial deposits are material of all grain sizes (clay, silt, sand, gravel, and cobbles) deposited by glacial meltwater streams in outwash plains and river valleys. Glaciolacustrine deposits generally consist of clay, silt, and fine sand deposited in temporary lakes that formed after the retreat of the glacial ice sheet. Stratified deposits are more widespread in eastern Massachusetts than in western Massachusetts. In eastern Massachusetts, stratified deposits include extensive outwash plains, particularly in the southeast (Stone and others, 2018). In other areas of the State, stratified deposits are more likely to be found in river valleys. On Cape Cod and the islands and areas of southeastern Massachusetts, the surficial geology is mainly stratified deposits (Stone and others, 2018). In these areas, precipitation percolates through the more permeable soils and unsaturated zone to the groundwater table (reducing surface runoff) and ultimately discharges to a pond or stream as base flow. Thus, runoff peaks in areas of extensive stratified deposit can be diminished in magnitude, and medium and lower flows generally have a higher component of base flow (in other words, groundwater discharge).

Hydrologic variability may also be associated with different physiographic provinces, and Denny (1982) identifies seven physiographic provinces within the study area. From eastern to western Massachusetts, the physiographic provinces are the Coastal Plain, coastal lowlands, central highlands, Connecticut Valley, Hudson-Green-Notre Dame highlands, Vermont Valley, and the Taconic highlands. Additionally, the EPA has divided the United States into ecological regions (U.S. Environmental Protection Agency, 2022b). These regions are based on ecosystems that generally are similar and have been identified through the analysis of the patterns and the composition of biotic and abiotic features. These features include geology, physiography, vegetation, climate, soils, land use, wildlife, and hydrology. The study area includes four EPA level III ecoregions: Atlantic Coastal Pine Barrens, Northeastern Coastal Zone, and Northeastern Highlands (U.S. Environmental Protection Agency, 2022b).

Land cover for the study area ranges from highly developed in and around cities in eastern Massachusetts, such as Boston (metropolitan area), Lowell-Lawrence, Brockton, Fall River, and New Bedford, to the less developed rural forested areas of communities in central and western parts of Massachusetts. However, central and western Massachusetts have highly developed areas in and around Worcester and Springfield, respectively, and several additional smaller cities. Overall, Massachusetts is about 64 percent forested and about 21 percent “built” (urban and suburban) (Harvard Forest, 2020). Water bodies and wetlands tend to be slightly more prevalent in eastern and southeastern Massachusetts, respectively, than in central and western Massachusetts—excluding Quabbin Reservoir in central Massachusetts.

Site Selection

All active and discontinued streamgages in Massachusetts with 8 or more years of record through September 30, 2022 (both water years¹ and climatic years² ), were evaluated for possible use in the regional regression analyses. Streamgages in Connecticut, Rhode Island, southern New Hampshire and Vermont, and eastern New York with at least part of their drainage areas within 25 miles of the Massachusetts border were also evaluated for the regression analyses. This list of streamgages included 174 streamgages (table 1); it excluded Mother Brook (01104000) (not shown) because it is a diversion channel and would not be used for at-site streamflow statistics or in regional regression equations for ungaged sites.

All potential streamgages were evaluated for flow regulations such as water withdrawals, diversions, flood control, hydropower generation, and wastewater discharge. Average annual withdrawal and wastewater discharge data for 2010–14 in Massachusetts (Levin and Granato, 2018) were retrieved from USGS StreamStats (https://streamstats.usgs.gov/ss/). Water-use data from Connecticut and Rhode Island were for annual water withdrawals and did not include wastewater discharge data (Laura Medalie, U.S. Geological Survey, written commun., 2021). No water withdrawal data were available for sites in New Hampshire, Vermont, and New York. Therefore, streamgages selected in these States were limited to those used in previous low-flow studies or known to have minimal to no regulations (Scott Olson and Andrew Waite, U.S. Geological Survey, oral commun., 2022). The EPA Enforcement and Compliance History Online (ECHO) database (U.S. Environmental Protection Agency, 2022a) was used for supplemental wastewater discharge data in the evaluations of the sites. Additionally, the USGS Gages II “hydrologic disturbance index” (Falcone, 2011) was used to evaluate sites on the basis of seven variables: (1) major dam density in 2009; (2) water withdrawals; (3) changes in dam storage, 1950–2009; (4) streams coded as a canal, ditch, or pipeline in the National Hydrography Dataset Plus (NHDPlus); (5) straight-line distance of the gage location to the nearest major NPDES point in the watershed; (6) road density; and (7) fragmentation index of undeveloped land in the watershed. Those streamgages with known regulations within their drainage basins that were substantial enough to clearly change the recorded daily mean streamflows due to dam operations, withdrawals, diversions, or wastewater discharges for more than several days during each water year were excluded from the dataset. Streamgages that have be used in previous low-flow studies were included in the site selection process. The final list of streamgages for the regression analyses included 81 streamgages with 10 or more climatic and water years of record and minimal to no regulations, located as follows:

Flow-Duration Statistics

Flow durations represent the percentage of time that a given flow is equaled or exceeded without regard to the sequence of recorded flows (Searcy, 1959). Typically, flow durations characterize the range of flow rates for the period over which data were collected. Flow durations were computed for complete water years for the entire period of record and for selected months for 174 streamgages (table 1 and fig. 1) with 8 or more complete water years of record in and near Massachusetts.

Flow durations are computed by sorting the daily mean streamflows for the period of interest (the entire record, a monthly period, or another period) from largest to smallest and assigning each streamflow value a rank, starting with one for the largest value. The frequencies of exceedance are then computed by using the Weibull plotting-position formula (Weibull, 1939):

Examples of flow-duration curves are provided for the streamgages Squannacook River near West Groton, Massachusetts (01096000; map number 12 in fig. 1 and table 1), and Mill River at Northampton, Mass. (01171500; map number 120 in fig. 1 and table 1), both of which are used in the mainland regression equations (fig. 3). Additional examples are provided for the streamgages Jones River at Kingston, Mass. (01105870; map number 58 in fig. 1 and table 1), and Quashnet River at Waquoit Village, Mass. (011058837; map number 61 in fig. 1 and table 1), both of which are used in the southeastern equations. Notably, Quashnet River at Waquoit Village, Mass. (011058837), has a very different shape to its flow-duration curves than the other streamgages. Its flow-duration curve is flatter (less of a slope), which is likely due to the surficial geology of its contributing drainage area. This streamgage’s drainage area is in the western part of Cape Cod and is mainly underlain by sand and gravel deposits (about 95 percent; Sturtevant and others, 2025). Its surface-water drainage area does not coincide with the groundwater contributing area, and the groundwater contributing area (10.6 mi²; Sturtevant and others, 2025) is significantly larger than the surface-water drainage area (2.58 mi²; Bent and others, 2025). Jones River at Kingston, Mass. (01105870), has a drainage area that is primarily underlain by sand and gravel deposits (about 77 percent, Bent and others, 2025) and is somewhat different than the other three flow-duration curves. This streamgage has a groundwater contributing area of 21.8 mi² (Sturtevant and others, 2025) and a surface-water drainage area of 20.1 mi² (both areas include the about 4.4 mi² contributing area to Silver Lake, a water supply for the City of Brockton, Bent and others, 2025). Differences in flow-duration curves can also be the results of different periods of record; regulations; and basin characteristics. Although these streamgages (fig. 3) have drainage areas ranging from about 2.58 to 65.9 mi² and their periods of records range from about 33 to 85 years for this study (table 1), the differences in their flow-duration curves are most likely due to the surficial geology.

The USGS Hydrologic Toolbox software was used to compute flow durations (Barlow and others, 2022). The selected flow durations range from 99 to 1 percent (table 2), with the number of selected durations increasing in the extreme percentile ranges (from 99 to 90 and from 10 to 1). Estimated flow statistics at the 174 streamgages in the study for their periods of record are available in Bent and others (2025).

Other streamflow statistics—annual flow durations (40-, 30-, 25, 20-, 15-, 10-, 5-, 2-, and 1-percent), monthly 90- and 50-percent flow durations (January–May and October–December), and median of the monthly means (January, March–May, July, and September–December)—were computed for streamgages for the mainland area (Bent and others, 2025) and for streamgages and partial-record stations in the southeastern area (Sturtevant and others, 2025).

Low-Flow Frequency Statistics

Low-flow frequencies are computed for streamgages by determining the frequency of an annual series for a consecutive number of low-flow days (Riggs, 1972)—for example, the 7-day, 10-year low flow-frequency (7Q10). This statistic is the minimum consecutive D-day mean streamflow that is expected to occur once in any Y-year period, or that has a probability of 1/Y of not being exceeded in any given year. Any combination of number of days of mean minimum streamflow and years of recurrence may be used to determine the low-flow frequencies. The annual series for the determination of low-flow frequencies for this study was based on a climatic year. Use of a climatic year rather than a water year allows for an analysis of an uninterrupted low-flow period; in Massachusetts, this low-flow period typically occurs from early August through mid-October. The minimum number of climatic years of record for the streamgages was 8 years, although all streamgages used in the mainland and southeastern regression equations had 10 or more climatic years.

For this study, low-flow frequencies were computed for the 7-day, 2-year (7Q2); 7-day, 10-year (7Q10); 30-day, 2-year (30Q2); and 30-day, 10-year (30Q10) statistics (Bent and others, 2025). Low-flow frequencies were computed by using the USGS Hydrologic Toolbox software (Barlow and others, 2022). An example plot of the annual 7-day low-flows with the log-Pearson type III distribution fit is shown in figure 4. The 7Q2 and 7Q10 streamflow are where, in figure 4, the annual non-exceedance probabilities of 50 and 10 percent, respectively, intersect with the log-Pearson type III curve.

Annual, Monthly, and Other Statistics

The streamflow statistics harmonic mean, monthly flow duration, median of the monthly means, and the median of the annual 7-day low flow were computed by using R packages (Bent and others, 2025). Harmonic mean was computed with the R DVstats package (U.S. Geological Survey, 2024b) according to the EPA DFLOW user’s manual (Rossman, 1990). The median of the annual 7-day low-flows was computed from the annual 7-day low-flow, which was determined by using the USGS Hydrologic Toolbox software (Barlow and others, 2022).

Trends in Low-Flows

The traditional assumption underlying low-flow analysis is stationarity in time. The assumption allows researchers to estimate low-flow statistics from past records and apply them to the future without adjustments. Milly and others (2008) called the assumption of climate-related stationarity into question and advocated for new methods to replace models based on stationarity. Several studies have shown that streamflows can be nonstationary by documenting increases in low and median flows across the United States (McCabe and Wolock, 2002; Lins and Slack, 2005; Small and others, 2006; Hodgkins and Dudley, 2011; Dudley and others, 2020).

For the trend analysis, the annual 7-day low-flow data were analyzed for long-term trends at streamgages in and near Massachusetts by using the same methodology as Ahearn and Hodgkins (2020). Subsets of streamgages with longer records were created to evaluate trends during the periods of 30, 50, 70, and 90 climatic years up to 2019 (tables 3, 4, 5, and 6, respectively). All 10-year blocks within each time period analyzed were required to be at least 80 percent complete so that no part of the time series would have substantial missing data. These length and completeness criteria resulted in 64 streamgages for the 30-year period (1990–2019), 58 streamgages for the 50-year period (1970–2019), 43 streamgages for the 70-year period (1950–2019), and 14 streamgages for the 90-year period (1930–2019). The magnitudes of the trends were computed with the Sen slope (also known as the Kendall-Theil robust line). The Sen slope was calculated by determining the median of all possible pairwise slopes in each time series (Helsel and Hirsch, 2002). The Sen slope is multiplied by the number of annual 7-day low flows to obtain the magnitude of the trend or total change in the annual 7-day low flows over the period analyzed. For example, a Sen slope of 0.099 cubic foot per second per year multiplied by 70 for the 70-year period results in a trend of 6.92 ft³/s for the North River at Shattuckville, Mass. (01169000; map number 114 in fig. 1 and table 1) (table 5).

The trends were computed with methods that consider the possibility of short-term persistence (STP) and long-term persistence (LTP) in the temporal data. This is an important issue that is often ignored in trend studies. Trends over time are sensitive to assumptions of whether underlying hydroclimatic data are independent, have STP, or have LTP (Cohn and Lins, 2005; Koutsoyiannis and Montanari, 2007; Hamed, 2008; Khaliq and others, 2009; Kumar and others, 2009). STP and LTP may represent the occurrence of wet or dry conditions that tend to cluster from year to year (Koutsoyiannis and Montanari, 2007; Hodgkins and others, 2017). For further discussion and references on persistence, refer to Hodgkins and Dudley (2011). Because the long-term time-series structure of low-flow data is not well understood, temporal trend significance with three different null hypotheses of the serial structure of the data is reported: independence, STP, and LTP (Hamed and Ramachandra Rao, 1998; Hamed, 2008). The serial structure of data referred to as “independence” means annual 7-day low flows from year to year are independent from each other (ignores any short or long clusters of wet and dry years). Trends were considered statistically significant at a p-value ≤0.05; this level represents a 5-percent probability that a trend is due to random chance. Results from the trend analyses for 30-, 50-, 70- and 90-year time periods under the three serial correlation structures, magnitudes of Sen slopes, and p-values are shown in tables 3 through 6. Low-flow trend results depend on the period of record analyzed and assumptions about the serial correlation structure of the annual peak flows.

For streamflow records influenced by regulation or other anthropogenic influences, interpretation of trend analyses is more complicated. Like near-natural sites, streamflow patterns at gages influenced by anthropogenic activities are also influenced by changes in climate patterns or basin characteristics. However, and this is especially true for regulated streamgages, those changes can be mitigated, enhanced, or even offset by changes in regulation patterns or other diversions. Nonetheless, trend assessments of flow patterns at such streamgages can be informative and, therefore, are included in these analyses.

For the 30-year period (1990–2019), 0 of the 24 streamgages used in the mainland regression analyses had statistically significant increasing or decreasing trends (p-value ≤0.05) if independence, STP, or LTP of 7-day annual low flows is assumed (table 3). For the southeastern regression analyses, one of two streamgages (Quashnet River at Waquoit Village, Mass.; 011058837; map number 61 in fig. 1 and table 1), had a statistically significant trend—increasing for independence and STP tests. For the other 38 streamgages not used in either of the regression analyses: 2 streamgages had a statistically significant increasing trend and 36 streamgages had no statistically significant trend for either independence, STP, or LTP tests.

For the 50-year period (1970–2019), 3 of the 21 streamgages used in the mainland regression analyses had statistically significant decreasing trends (p-value ≤0.05) (Nashoba Brook near Acton, Mass., 01097300, map number 17; Parker River at Byfield, Mass., 01101000, map number 26; and Branch River at Forestdale, Rhode Island, 01111500, map number 80—in fig. 1 and table 1) if independence, STP, or LTP of 7-day annual low flows is assumed (table 4). For the one streamgage in the southeastern regression analyses, there was no statistically significant trend. For the other 36 streamgages not used in either of the regression analyses: 3 streamgages had a statistically significant increasing trend, 5 had a decreasing trend, and 28 had no statistically significant trend in either independence, STP, or LTP tests. Two of these streamgages not used in either regression analyses also had statistically significant increasing trends for the 30-year period: Swift River at West Ware, Mass. (01175500; map number 134 in fig. 1 and table 1), and Deerfield River at Charlemont, Mass. (01168500; map number 113 in fig. 1 and table 1). Swift River at West Ware, Mass. (01175500), is downstream from Quabbin Reservoir, which is the water-supply for much of the Boston metropolitan area, and flows at Deerfield River at Charlemont (01168500) are affected by hydropower generation.

For the 70-year period (1950–2019), 4 of the 13 streamgages used in the mainland regression analyses had statically significant trends (p-value ≤0.05): 3 streamgages increasing and 1 streamgage decreasing if independence, STP, or LTP of 7-day annual low flows is assumed (table 5). No streamgages used in the southeastern regression analyses have period of records that go back 70 years. For the other 30 streamgages not used in either of the regression analyses: 6 streamgages had a statistically significant increasing trend, 2 had a decreasing trend, and 22 streamgages had no statistically significant trend in either independence, STP, or LTP tests. Only one of the streamgages also had statistically significant similar trends for the 30- and 50-year periods: increasing for Deerfield River at Charlemont, Mass. (01168500; map number 113 in fig. 1 and table 1).

For the 90-year period (1930–2019), only one streamgage was used in the mainland regression analyses, and it did not have a statistically significant trend for the other three assessments (table 6). No streamgages used in the southeastern regression analyses have periods of records that go back 90 years. For the other 13 streamgages not used in either regression analyses, 2 had statistically significant trends: 2 increasing and 2 decreasing in either independence, STP, or LTP tests. Only one of these streamgages also had similar statically significant trend for the 30-, 50-, and 70-year periods: increasing for Deerfield River at Charlemont, Mass. (01168500; map number 113 in fig. 1 and table 1), where flows are affected by hydropower generation.

As the science evolves and new data are obtained, further analysis could improve understanding of the trends observed in this study and the effects on low flows. Historical low-flow trends in and near Massachusetts do not offer clear and convincing evidence of the need to incorporate trends into low-flow analyses. If the evidence becomes clear, a well-defined deterministic mechanism should be identified prior to incorporating trends (Salas and others, 2018). For this study, the traditional assumption of stationarity is used with no adjustment for historical trends.

Basin Characteristics

The characteristics of streamflow are directly related to a drainage basin’s physical, land-cover, land-use, geologic, and climatological characteristics (table 7). Characteristics of the drainage basin were selected for use as potential explanatory variables in the regression analysis on the basis of their theoretical relations to low flows, the results of previous low-flow studies in similar hydrologic regions, and the feasibility of determining the basin characteristics with digital datasets and GIS technology. Measuring the basin characteristics with GIS technology facilitates automation of the process of solving the regional regression equations by using the USGS StreamStats web-based application.

The basin boundaries delineated with StreamStats for the 174 streamgages in and near Massachusetts were overlaid on areal coverages of the basin characteristics of interest to determine the characteristics’ values for the basin upstream from each site (Bent and others, 2025). Basin, land-use, land cover, surficial geology, soil, and climatological characteristics were determined for the 174 streamgages in and near Massachusetts. These data and the sources of the GIS data are published in a USGS data release (Bent and others, 2025).

Streamflow Variability Index

Streamflow variability index (SVI) is a measure of the variability in streamflow resulting from variability in precipitation, as mitigated by characteristics of the basin such as surface storage and groundwater discharge (base flow). SVI has been found to be an explanatory variable in low-flow equations for recent studies in Alabama (Feaster and others, 2020), Iowa (Eash and Barnes, 2012), Kentucky (Martin and Ruhl, 1993; Martin and Arihood, 2010), Missouri (Southard, 2013), Ohio (Koltun and Whitehead, 2002; Whitehead, 2002; Koltun and Kula, 2013; VonIns and Koltun, 2024), and West Virginia (Friel and others, 1989). Unregulated streams with relatively small SVIs tend to have proportionally more flow contributed from groundwater discharge and (or) surface storage than streams with larger SVIs (Searcy, 1959). Figure 3 shows the flow duration curves for four selected streamgages and their associated calculated SVIs. Note that the Quashnet River at Waquoit Village, Mass. (011058837; map number 61 in fig. 1 and table 1), has a small SVI and has a relatively high more contribution to flow from groundwater discharge given that 95 percent of the drainage area is underlain by sand and gravel deposits (Sturtevant and others, 2025). Lane and Lei (1950) proposed SVI as a method to help produce synthetic flow-duration curves. The SVI is defined as the standard deviation of the logarithms of the 19 streamflow values at 5-percent class intervals from 5 to 95 percent on the daily flow-duration curve for the analysis period (Searcy, 1959). The formula for the SVI discussed in this report is

SVIs determined initially from streamgages with 8 or more water years of record in southern New England and eastern New York were plotted on a map (not shown) to assess spatial trends. Although there were visually identifiable spatial trends (for example, a cluster of low SVIs at streamgages in southeastern Massachusetts and Cape Cod—an area known for relatively high groundwater discharge)—it was apparent that, in some areas, SVIs can change appreciably over relatively small distances of 10–20 miles. Consequently, it was deemed important to compute and use as much SVI data as possible to prepare the grid. Therefore, in development of an SVI map for southern New England and eastern New York, SVI was computed at additional streamgages (some with periods of record less than 8 water years) and partial-record stations to improve the SVI map (Bent and others, 2025).

Koltun and Kula (2013) also estimated SVIs for other streamgages and partial-record stations in Ohio to assist in development of a detailed SVI map. Streamgages and partial-record stations within southern New England and eastern New York with published flow durations were added to the SVI database for creating the map (Bent 1995, table 5; Ries 1999, table 3; Bent, 1999, tables 8 and 9; Bent and others, 2014, tables 3 and 6). For the streamgages with a period of record less than 8 years, the flow-duration curve was used to compute the SVI for that streamgage. But for most of the partial-record stations, only flow durations from the 99th to 50th percentiles were available because they were mainly low-flow partial-record stations. Therefore, a relation between streamgages’ SVIs for the 10 flow durations in 5-percent class intervals from the 95th to 50th percentiles (SVI_95–50) and the 19 flow durations in 5-percent class intervals from the 95th to 5th percentiles (SVI_95–5) was developed (fig. 5). Development of this relation started with 304 streamgages with 8 or more water years of record in southern New England and eastern New York. During analyses, 56 streamgages were removed from the relation as they were determined to be outliers. Evaluation of these 56 streamgages’ flow-duration curves from the 99th to 1st percentiles showed clear changes that were associated with regulations, such as major water supply withdrawals, flood-control dams, and hydropower generation. The final relation between the SVI_95–50 and the SVI_95–5 was determined by using 248 long-term streamgages in southern New England and eastern New York (Bent and others, 2025). The relation is very good, with a coefficient of determination (R²) of 0.9583. The power regression equation for this relation was then used to estimate the SVI for the 95th- to 5th-percentile flow durations at partial-record stations that only have published flow durations from the 95th to 50th percentiles in 5-percent class intervals:

SVI must be regionalized in some fashion to permit its estimation at ungaged sites. The geostatistical techniques of kriging, inverse distance weighting, and natural neighbor were all examined as possible methods for regionalizing SVI (Esri, 2023c). Maps of SVI regionalized by geostatistical techniques were developed by using 664 sites (276 streamgages with 8 or more water years of record, 102 streamgages with less than 8 water years of record, and 286 partial-record stations) in and near southern New England and eastern New York (Bent and others, 2025). Of the 276 streamgages with 8 or more water years of record, 28 streamgages were just outside the southern New England and eastern New York area. These streamgages were included to make sure that the drainage areas of streamgages on the perimeter of the area being used in the regression analyses would be covered by the map to determine a mean basin SVI. Evaluation of maps of the regional SVI for southern New England and eastern New York created using the geostatistical techniques determined kriging to be the most representative method (Esri, 2023b). Kriging is based on statistical models that include autocorrelation, which is the statistical relations among the measured data points (Esri, 2023a). Kriging fits a specified number of data points within a specified radius to determine output values. It has multiple steps, including exploratory statistical analyses, variogram modeling, and creating the output values. Kriging is often used when there is spatial correlation or bias in the data, especially with soil and geologic related data. Thus, kriging is appropriate for SVI data. The final SVI grid for southern New England and eastern New York is illustrated in figure 6. This SVI grid was used to determine a mean SVI value for each streamgage’s drainage area in and near southern New England. The mean basin SVI values for the 174 streamgages tested for possible use in the regression analyses are in Bent and others (2025).

Comparison of the mean SVI values from kriging to the computed at-site SVI values using equation 2 was possible at 135 of the 174 streamgages. Of the 39 streamgages not used in the comparison, 37 sites had regulations, such as major water supply withdrawals, flood-control dams, and hydropower generation; and for 2 sites on Cape Cod, an accurate surface-water drainage area could not be determined. The average and median of the mean SVI values were 0.468 and 0.464, respectively, which compared well to the average and median of the at-site SVI values of 0.472 and 0.464, respectively, for the 135 streamgages.

Drainage-Area Ratio Methods

The drainage-area ratio method is based on the assumption that the near-natural (minimal to no regulations) streamflow at an ungaged site is the same per unit of drainage-basin area as that at a nearby hydrologically similar streamgage with near-natural flows. Drainage areas of the ungaged site and the streamgage are determined from topographic maps or the USGS StreamStats web-based application (https://streamstats.usgs.gov/ss/). Streamflow statistics are computed for the streamgage, and then the statistics (streamflow values) are divided by the drainage area to determine the streamflow for each statistic per unit area, in cubic feet per second per square mile, at the streamgage. These values are multiplied by the drainage area to the ungaged site to obtain estimated statistics for that site. This method is most commonly applied if the index streamgage is on the same stream as the ungaged site because the accuracy of the method depends on the proximity of the two sites, on similarities in drainage area, and on other physical and climatological characteristics of their drainage basins.

Several studies have provided estimates of the maximum difference in drainage areas for which the use of the drainage-area ratio method would generate more accurate estimates of streamflow statistics than the use of regression equations. Guidelines have been provided for estimating peak-flow statistics, and usually the recommendation has been that the drainage area to the ungaged site should be 0.5 to 1.5 times the drainage area of the index streamgage (Choquette, 1988, p. 41; Koltun and Roberts, 1990, p. 6; Lumia, 1991, p. 34; Bisese, 1995, p. 13; Koltun and Whitehead, 2002, p. 22; Martin and Arihood, 2010, p. 28). Koltun and Schwartz (1987, p. 32) recommended a narrower range of 0.85 to 1.15 times the drainage area of the index streamgage for estimating low flows at ungaged sites in Ohio. Ries and Friesz (2000), however, determined that the drainage-area ratio method could be used to estimate low-flow statistics for ungaged sites in Massachusetts if the drainage area for an ungaged site was between 0.3 and 1.5 times the drainage area of the index streamgage site. They found that this method was generally as accurate as, or more accurate than, regression equations for this range of drainage areas.

In the drainage-area ratio method, the streamflow values are transferred from a streamgage to the ungaged site by the following formula:

Eash and Barnes (2012) compared estimates of the 7-day, 10-year low flow developed by using regional regression equations, the drainage-area ratio method, and the weighted drainage-area ratio method to estimates from streamflow records for 48 streamgages (31 pairs of streamgages) on rivers in Iowa. They found that the weighted drainage-area ratio provided the best estimate if the drainage area to the ungaged site was between 0.4 and 1.5 times the drainage area of the streamgage.

In the weighted drainage-area ratio method, the streamflow values are transferred from a streamgage to the ungaged site by using the following formula:

As the ratio (Q_gs/Q_gr) approaches 1, or the ratio of DRNAREA_u to DRNAREA_g approaches 0.5 or 1.5, the weighting factor in equation 5 approaches 1 and it no longer has an effect on the regression equation estimate for the ungaged site (Q_ur). Additionally, both the drainage-area and weighted drainage-area ratio methods may not be applicable for ungaged sites where physical, land-cover, land-use, surficial-deposit, or climatological characteristics or regulations are substantially different between the ungaged site and the streamgage. The error associated with estimates based on the drainage-area and weighted drainage-area ratio methods cannot be calculated.

Regional Regression Analysis

Studies to develop regional regression equations for estimating streamflow statistics at ungaged sites have been done in many States throughout the United States, including those adjacent to Massachusetts (refer to the report section “Previous Studies”). Multiple regression analyses provide a mathematical equation for estimating a response (dependent) variable—that is, a streamflow statistic—from one or more explanatory (independent) variable(s), such as basin characteristics. Ideally, the development of regression equations involves the use of streamflow data from a large number of long-term streamgages on unaltered streams evenly distributed across the region of interest and with a range of basin characteristics. But in many cases, the number of streamgages representing unaltered flow is limited, the streamgage network is biased toward representing larger streams or rivers, the network is unevenly distributed geographically, and the range of basin characteristics upstream from streamgages does not cover the complete range found in the region. In Massachusetts, the number of streamgages on streams with minimal to no regulations is limited; for this reason, nearby streamgages in the surrounding States within about 25 miles of the border were used. The daily mean streamflow records at these streamgages were used to compute selected low-flow statistics, and then those statistics and basin characteristics at these streamgages were used to develop regional regression equations for estimating selected low-flow statistics at ungaged sites in Massachusetts.

Multiple regression is used to create equations that relate streamflow statistics for streamgages to the physical, land-cover, land-use, surficial-deposit, soil, and climatological characteristics of their upstream drainage areas. Once an optimal equation has been determined, a streamflow statistic at a nearby ungaged site in a basin with similar characteristics can then be estimated by applying the equation to the ungaged site.

The basic equation describing a linear multiple regression analysis is

Linear regression analysis is based on the following assumptions: (1) the mean of the residuals (e_i) is zero, (2) the variance of the residuals is constant, (3) the residuals are normally distributed, and (4) the residuals are independent of each other. In addition to these assumptions, the selected explanatory variables (X) should have a physical basis as predictors of the streamflow statistic, the explanatory variables (basin characteristics) in the equation should not be highly correlated with each other, and the signs of the terms of the equation should make hydrological sense. For example, the variable drainage area should have a positive coefficient because an increase in drainage area should result in an increase in the value of the streamflow statistic.

In almost all regionalization studies to determine low-flow and peak-flow statistics, the response and explanatory variable datasets are skewed. As a consequence, the data need to be transformed to ensure that the mean of the residuals equals zero. In many studies, a logarithmic transformation is used. A base-10 log-transformed multiple regression equation has the form

After the coefficients have been determined through regression analysis, the equation is transformed back to its original units in a form that can be used to estimate a specific streamflow statistic at an ungaged site. The retransformed equation has the following form:

A linear regression equation provides an unbiased estimate of the mean response of the response variable. Although estimates provided by equation 7 are unbiased, these estimates are in log units, whereas estimates in the original units are needed to calculate specific streamflow statistics at an ungaged site. Estimates from equation 8 are in the original units. However, this equation predicts the median, instead of the mean, value of the response variable. A streamflow statistic based on a median creates an estimate that is biased and tends to be lower than the mean (Ries and Friesz, 2000). Bias correction factors (BCFs) were used in some studies in Massachusetts and New Hampshire to remove the bias from the estimate (Ries, 1994a, b; Ries and Friesz, 2000; Flynn, 2003a, b; Archfield and others, 2010). In other studies (Risley, 1994; Stuckey, 2006; Armstrong and others, 2008; Ahearn, 2010), BCFs were not used because they were generally very small. In this study, BCFs were also not used because if they were, then the streamflows estimated from the regression equations would not have an equal chance of being higher or lower than their actual values (Julie Kiang, U.S. Geological Survey, oral commun., 2011).

Final Regression Equations

The final regression equations were determined by using the technique in the weighted-multiple-linear regression (WREG) program (Eng and others, 2009) and left-censored regression (Kroll and Stedinger, 1999). This was done by using the R package for WREG (Farmer, 2023) and the R package for left-censored regression (Lorenz, 2015).

The final regression equations for the flow-duration statistics from 99 to 50 percent; June to September 90- and 50-percent flow duration; median of the February, June, and August means; harmonic mean; and median of the annual 7-day low flow were developed by using weighted least-squares (WLS) regression. WLS is used for these low-flow statistics because they are not frequency statistics that involve the log-Pearson type III distribution for their estimation.

For the low-flow frequency statistics 30Q2 and 30Q10, the final regression equations were developed by using generalized least-squares regression (GLS). GLS requires the skew, K, and standard deviation for the log-Pearson type III distribution (Eng and others, 2009), which were computed by using the USGS Hydrologic Toolbox (Barlow and others, 2022).

For the low-flow frequency statistics 7Q2 and 7Q10, the final regression equations were developed by using weighted left-censored regression. This was because, for the 7Q2 statistic, 1 of 81 streamgages had a flow value of zero and, for the 7Q2 statistic, 5 streamgages had flow values of zero. The value of zero cannot be logarithmically transformed, but logarithmic transformation of streamflow statistics (that are not normally distributed) is needed to develop regression equations. The method of handling the zero values in a regression analysis depends on the number of streamgages in the dataset with response variables equal to zero. If the number is small, left-censored regression analysis is used: either an adjusted maximum-likelihood estimation (Cohn, 1988; Helsel, 2005) or a censoring method that adds a small constant value (0.01, 0.10, and 1.00) to all response variables in a dataset (Kroll and Stedinger, 1999). Censoring and coding data as “less than” a threshold value (0.01, 0.10, and 1.00) allows the use of a log transformation on the data and, therefore, allows all the data (uncensored and censored) to be used in the regression analysis to develop the regional equations (Watson and McHugh, 2014). For the 7Q2 and 7Q10 regression analyses, a small constant value of 0.01 (censoring threshold value of 0.01 ft³/s) was used because this is the lowest flow value reported by the USGS. Weighted left-censored regression techniques were used for the final regression analyses, with the weights being based on the number of years of record at the streamgages (Ziegeweid and others, 2015; Gotvald, 2017; Feaster and others, 2020).

The regression analyses determined that drainage area (represented by the variable DRNAREA) was a significant explanatory variable (p-value less than or equal to 0.05) in all 27 of the regression equations. Drainage area generally is the most significant explanatory variable in all regional streamflow regression equations, whether for low flows; peak flows; or mean annual, mean monthly, or median monthly statistics. The final regression equations for the 27 streamflow statistics are listed in table 8.

Table 8.    

Summary of regression equations and measures of model accuracy for estimating selected streamflow statistics for the mainland area of Massachusetts (excludes the Plymouth-Carver-Kingston-Duxbury aquifer system area in southeastern Massachusetts and Cape Cod).

[Flow duration is the percentage of time that streamflow was equaled or exceeded. R², coefficient of determination; MSE, mean square error in log-base 10 logarithm; RMSE, root mean square error, in percent; DRNAREA, drainage area in square miles; SOILAB+0.1, percent area of the combined Soil Survey Geographic Database (SSURGO; Natural Resources Conservation Service, 2022) hydrologic soil type A and B plus the value of 0.1 to ensure there are no zero values; WLS, weighted least-squares regression; WREG, weighted-multiple-linear regression (Eng and others, 2009); SVI, streamflow variability index (unitless); TEMP, mean annual temperature; GLS, generalized least-squares regression]

Table 8.    Summary of regression equations and measures of model accuracy for estimating selected streamflow statistics for the mainland area of Massachusetts (excludes the Plymouth-Carver-Kingston-Duxbury aquifer system area in southeastern Massachusetts and Cape Cod).

Streamflow statistic	Regression equation	Regression method	Number of streamgages	Adjusted R2	MSE, logarithm base-10	RMSE, in percent
Flow duration
50	0.67788(DRNAREA)1.0090(SOILAB+0.1)0.1306	WLS–WREG	81	0.9846	0.0059	17.85
60	0.46212(DRNAREA)1.0248(SOILAB+0.1)0.1345	WLS–WREG	81	0.9808	0.0077	20.35
70	0.26802(DRNAREA)1.0612(SOILAB+0.1)0.1424	WLS–WREG	81	0.9716	0.0124	26.11
75	0.10816(DRNAREA)1.0668(SVI)−1.4915	WLS–WREG	81	0.9616	0.0181	31.70
80	0.05675(DRNAREA)1.0957(SVI)−1.8698	WLS–WREG	81	0.9505	0.0252	37.78
85	0.02824(DRNAREA)1.1285(SVI)−2.2595	WLS–WREG	81	0.9346	0.0363	46.03
90	0.01215(DRNAREA)1.1685(SVI)−2.7577	WLS–WREG	81	0.9153	0.0518	56.20
95	0.00381(DRNAREA)1.2276(SVI)−3.4435	WLS–WREG	81	0.8809	0.0850	75.45
98	0.00113(DRNAREA)1.2923(SVI)−4.1659	WLS–WREG	81	0.8380	0.1370	103.33
99	0.000458(DRNAREA)1.3474(SVI)−4.7437	WLS–WREG	81	0.8068	0.1879	130.70
Monthly flow duration
June 50	0.38485(DRNAREA)1.0436(SOILAB+0.1)0.1750	WLS–WREG	81	0.9741	0.0106	24.08
July 50	0.05342(DRNAREA)1.1022(SVI)−2.0592	WLS–WREG	81	0.9430	0.0295	41.17
August 50	0.02557(DRNAREA)1.1399(SVI)−2.3736	WLS–WREG	81	0.9197	0.0468	53.06
September 50	0.01980(DRNAREA)1.1367(SVI)−2.6875	WLS–WREG	81	0.9261	0.0425	50.27
June 90	0.08364(DRNAREA)1.1291(SOILAB+0.1)0.2419	WLS–WREG	81	0.9319	0.0346	44.85
July 90	0.00590(DRNAREA)1.2300(SVI)−3.0001	WLS–WREG	81	0.8839	0.0795	72.42
August 90	0.00126(DRNAREA)1.3049(SVI)−4.0821	WLS–WREG	81	0.8381	0.1406	105.22
September 90	0.001128(DRNAREA)1.2786(SVI)−4.2534	WLS–WREG	81	0.8293	0.1413	105.58
Median of the monthly means
February	0.000000013(DRNAREA)1.0085(SOILAB+0.1)0.0689(TEMP)4.8020	User WLS–WREG	81	0.9813	0.0070	19.38
June	0.58275(DRNAREA)1.0336(SOILAB+0.1)0.1279	User WLS–WREG	81	0.9762	0.0096	22.89
August	0.07521(DRNAREA)1.0689(SVI)−1.6675	User WLS–WREG	81	0.9202	0.0408	49.11
Frequency
7Q2	0.00388(DRNAREA)1.2640(SVI)−3.2570	Weighted left-Censored1	81	20.9089	0.2601	65.69
7Q10	0.000310(DRNAREA)1.3850(SVI)−4.8170	Weighted left-Censored1	81	20.8314	0.4143	121.80
30Q2	0.019720(DRNAREA)1.1680(SVI)−2.1050	GLS–WREG	81	0.9093	0.0558	58.68
30Q10	0.001786(DRNAREA)1.3180(SVI)−3.4250	GLS–WREG	81	0.8217	0.1673	119.47
Other
Harmonic mean	0.02311(DRNAREA)1.1761(SVI)−2.9982	WLS–WREG	81	0.9003	0.0622	62.53
Median of the annual 7-day low-flow	0.00353(DRNAREA)1.2659(SVI)−3.3843	WLS–WREG	81	0.8663	0.1050	86.32

¹

Left-censored regression using a threshold of 0.01 cubic feet per second.

²

Psuedo R².

The SVI (unitless) was a significant explanatory variable (p-value less than or equal to 0.05) in 20 of the 27 regression equations and was in the equations estimating lower streamflow statistics (for example, the 99- to 75-percent flow durations). The combined percentage of SSURGO hydrologic soil groups A and B (SOILAB) was a significant explanatory variable in the remaining 7 of the 27 regression equations, and this variable tended to be in the equations estimating more moderate streamflow statistics (for example, the 70- to 50-percent flow durations). Mean annual temperature from 1981 to 2010 (TEMP) was in one equation, to estimate the February median of the monthly means.

As noted previously, several recent low-flow studies have SVI as an explanatory variable for low-flow statistic regression equations—for example, in Alabama, Iowa, Kentucky, Ohio, and West Virginia. Both hydrologic soil groups A and B were explanatory variables in low-flow equations in Iowa (Eash and Barnes, 2012).

The coefficient for SVI was negative in the 20 regression equations it was in, and the coefficient was more negative for the lower streamflow statistics, −4.7437 coefficient at the 99-percent flow duration and −1.4915 coefficient at the 75-percent flow duration (table 8). As shown in figure 3, streamgages with larger SVIs (for example, Squannacook River near West Groton, Mass.; 01096000; map number 12 in fig. 1 and table 1) tend to have a flow-duration curve with a steeper slope. In order to reduce the larger SVI values, the coefficient needs to become more negative as the streamflows get lower.

The coefficient SOILAB was positive in all the regression equations, which is expected because hydrologic soil group A consists of well-drained sands and gravelly sands with high infiltration and low runoff rates and group B consists of well-drained, moderately fine to moderately coarse-textured soils with moderate rates of infiltration and runoff (Esri, 2022). SOILAB tends to be located in areas with more sand and gravel, alluvium, and fine sand surficial deposits, and therefore those areas tend to have higher contributions from groundwater discharge (base flow) at the lower streamflows. SOILAB is shown in the regression equations as SOILAB+0.1 (table 8). This is to ensure that no SOILAB value could be zero, because a value of zero cannot be logarithmically transformed.

Mean annual temperature (TEMP) was a positive coefficient in the regression equation for the February median of the monthly means (table 8). During February, those areas with higher mean annual temperatures may have more runoff, due either to more precipitation as rain rather than snow or to more snowmelt, than areas with slightly colder mean annual temperatures.

The RMSE for the 27 regression equations ranged from 17.85 to 130.70 percent (table 8). The low-flow statistics—the 99- and 98-percent flow durations, August and September 90-percent flow durations, and 7Q10 and 30Q10—had standard errors of the estimate greater than 100 percent. Generally, the higher the flow statistic, the lower the standard error of the estimate. Although the standard errors of the estimate are high for the lower flow statistics, they were similar to those calculated for regional regression equations developed to estimate low-flow statistics in Rhode Island (Bent and others, 2014), Connecticut (Ahearn, 2010), and Massachusetts (Ries and Friesz, 2000), and in other USGS low-flow studies.

Review of areal plots (not shown) of the residuals (differences between streamflow statistics estimated from measured streamflow and those estimated from the regression equations) for selected low-flow frequencies and durations at the 81 streamgages did not indicate any strong regional biases (clear groups of negative or positive residuals). Additionally, plots of the streamflow statistics estimated from measured streamflow (observed data) as functions of the same statistics estimated from the regression equations (predicted data) are presented in figures 7A–G. These plots show no clear bias of the streamflow values generated by regression equations with respect to measured streamflows versus predicted values. However, at lower flows (7Q10, 99-percent flow duration, 30Q2, and 90-percent flow duration), the range of the data is clearly wider than that for slightly higher low-flows (median of the monthly means for August, harmonic mean, and June 50-percent flow duration). Again, a wide range between observed and predicted data is generally found in all low-flow studies, especially for the equations for estimating lower flows with the higher percent RMSE—for example, equations where the RMSE is greater than 50 percent (table 8).