Stationary Flood Frequency

Flood frequency analysis uses statistical techniques to estimate flood magnitudes associated with specific AEPs. An AEP is the probability that a flood of a specific magnitude or higher will occur in a given year. Flood frequency estimates at a streamgage are computed by fitting a probability distribution to the time series of the logarithm of annual peak streamflows. Flood magnitudes for AEPs of interest are computed from the quantiles of the fitted probability distribution. Bulletin 17C recommends fitting a Log-Pearson type III (LP3) distribution for peak streamflow (England and others, 2018). For this study, the method of moments was used to estimate the parameters of the LP3 distribution and corresponding quantiles. This method uses the sample mean (m_y), standard deviation (s_y), and skew (g_y) of the logarithms of the peak streamflows (y) to estimate distribution moments, which is shown in equations 1–3. The lmomco package in R (Asquith, 2023) was used to fit the LP3 distributions and compute quantiles from the sample moments.

The quantity (1+6/n) is a bias correction factor for the skew of an LP3 distribution recommended by Tasker and Stedinger (1986). The stationary FFA method used in this study differs from that recommended by Bulletin 17C, which uses the expected moments algorithm to fit the distribution and also recommends using a weighted or regional skew (England and others, 2018). The method of moments using the sample skew is used here so that comparisons between stationary and nonstationary estimates of flood magnitudes can be made with more consistent methodologies. Bulletin 17C recommends using the multiple Grubbs-Beck test (Cohn and others, 2013) to identify and potentially remove low outlier values to improve the fit of the distribution for the highest streamflows. This method, which assumes stationarity, may not be suitable in NFFA because the thresholds for identifying an outlier may change as the peak streamflow distribution changes through time. Because there is currently no recommended method for identifying or handling outliers in NFFA, and to keep stationary and nonstationary FFA methodologies consistent, all outliers were retained in the dataset.

Conditional Peak Streamflow Moments

When peak streamflows exhibit nonstationarity, the peak streamflow probability distribution changes over time. To estimate a flood AEP for current conditions, the conditional peak streamflow distribution for that year of interest must be determined. Serago and Vogel (2018) introduce expressions for estimating conditional distribution moments derived from a linear regression equation relating nonstationary peak streamflows to an explanatory variable such as a time index, urbanization, or climate variable (eq. 4).

If the regression residuals are homoscedastic and normally distributed, the conditional mean of the peak streamflow distribution can be obtained by solving the regression equation for a specific value of ω (eq. 5).

The conditional standard deviation can be derived by substituting equation 5 into formulas for standard deviation and skew (eqs. 2 and 3).

The difference $y_i-{\mu }_{y|{\omega }_i}$ is equivalent to the residuals of the fitted regression model and p adjusts for the degrees of freedom in the regression if more than one explanatory variable is used.

There are several derivations for conditional skew based on an OLS regression. Serago and Vogel (2018) and Hecht and others (2022) derive conditional skews based on an OLS regression using time as an explanatory variable. The simplifications made in these derivations apply only to uniformly distributed variables and are not applicable to a stochastic climate variable. This report uses an alternative method of computing the skew introduced in Glas and others (2023), which is applicable for time and climate explanatory variables. This method computes the skew using a standardized time series of peak streamflows:

The conditional skew (γ_Z) can then be computed using Z_i instead of y_i in equation 3:

Equations 6 and 8 produce a conditional standard deviation and skew that are constants and therefore are most appropriate when the nonstationarity arises from a change in the mean only.

Regression equations used for computing conditional peak streamflow moments should have homoscedastic and normally distributed residuals to avoid biased estimates of conditional moments. Regression residual assumptions were tested using a Breusch-Pagan test of homoscedasticity (Breusch and Pagan, 1979) and a probability plot correlation coefficient (ppcc) test for residual normality (Vogel, 1986). The Breusch-Pagan test was computed using the bptest function in the lmtest package in R (Zeileis and Hothorn, 2002). The null hypothesis for this test is that residuals are homoscedastic. Probability values (p-values) greater than 0.05 indicate an acceptable level of homoscedasticity for regression models. The ppcc test was computed in R using the ppccTest function in the ppcc package (Pohlert, 2020a). This test compares the correlation between the ordered residuals and the theoretical quantiles of a normal distribution. P-values greater than 0.05 and correlations close to 1 indicate that the residuals are likely to be normally distributed.

If the mean and variability of peak streamflows are nonstationary, the regression equation used to model the conditional distribution moments may have heteroscedastic residuals. Hecht and Vogel (2020) developed an approach to estimate the conditional standard deviation (${\sigma }_{y|{\omega }_0})$for a regression equation with heteroscedastic residuals. This approach uses a second regression to statistically model the variance of the fitted regression equation (eq. 9), which is then used to determine the conditional variance of the peak streamflow distribution (eq. 10). The residuals, ε, from equation 4 are used as the dependent variable in the variance regression. Equation 9 relates an independent variable to the residuals that have been raised to the 2/3 power. This transformation, called an Anscombe transformation (Anscombe, 1961), produces normally distributed, nonnegative values.

If the residuals of equation 9 are normally distributed and homoscedastic, the conditional variance of peak streamflows (y) given ω₀ can be computed according to equation 10 (refer to Hecht and Vogel [2020] for a complete derivation).

The conditional standard deviation, ${\sigma }_{y|{\omega }_0}$, in the heteroscedastic case is the square root of equation 10. For this study, the explanatory variable used in equation 10 was assumed to the same as that used in equation 4, although this is not a requirement of the method. If there is evidence that there are different causal mechanisms for the change in mean and the change in variance of peak streamflows, they could be modeled using different explanatory variables (Hecht and Vogel, 2020).

The conditional moments (eqs. 5, 6, 8, 10) can be used to fit a conditional peak flow distribution using the method of moments in the same way as for stationary FFA. Conditional moments can be used to fit any of the common distributions used in FFA. In this study the LP3 distribution is used to maintain consistency with the stationary estimates.

Flood Probability Conditioned on a Stochastic Variable

When the explanatory variable used to model peak streamflows is time, the flood magnitude for a given AEP can be computed directly from the conditional distribution. In this case, the conditional AEP represents the probability that a flood of a given value will occur in year ω₀. For example, to update the 1-percent AEP flood estimate to conditions in the year 2020 using a regression equation with water year as the explanatory variable (ω), the conditional mean, standard deviation, and skew could be computed for ω₀=2020. Then, the 1-percent AEP could be determined from the quantile of the fitted distribution.

If a climate variable is used as the explanatory variable, ω, in the regression equation (eq. 4), the conditional peak streamflow moments can be computed for a specific value of ω₀ using equations 5–10. To estimate the flood frequency for conditions in a specific year of interest, an appropriate value of ω₀ is needed for the year of interest. However, because climate variables have random variability from year to year, there is a range of plausible values that ω₀ could take in any given year.

For example, annual precipitation and peak streamflow at USGS streamgage 05481000 (Boone River near Webster City, Iowa) have upward trends in peak streamflow and an upward trend in precipitation (fig. 3A, peak streamflow not shown). A conditional flood frequency curve can be generated for peak streamflow for any single value of precipitation using the relation between peak streamflow and precipitation (fig. 3B). When estimating the flood frequency for the year 2020 using equations 5–10, a representative value of precipitation (ω₀) is needed. In 2020, the observed annual precipitation was 27.8 inches, the mean annual precipitation estimated by a linear regression is 34.0 inches, and the median annual precipitation estimated by a quantile regression is 35.8 inches. Any of these values may be a plausible value to use for ω₀ in equation 5; however, the conditional mean streamflow and the resulting AEP flood magnitude estimate will vary considerably depending on the choice of ω₀ (fig. 3A, B).

The conditional probability of a flood for a specific value of annual precipitation is of limited use in most applications because there is a range of possible annual precipitation values that could occur in any given year. Instead, the conditional peak streamflow distribution can be weighted by the probability distribution of annual precipitation for the year of interest to produce a probability-weighted estimate of the flood frequency for the year 2020 (fig. 3C).

The total probability, P(y), of a flood across all probable values of a climate variable $\omega$ can be computed using the law of total probability (Harchol-Balter, 2024).

In terms of the NFFA, this law can be generally interpreted as a weighted average of the conditional probability of a flood magnitude (y) given some value of precipitation, weighted by the relative probability of that annual precipitation value, and evaluated across the whole distribution of annual precipitation values.

When obtaining the total probability of a flood of a given magnitude using equation 11, it is necessary to know the probability density of the explanatory variable used to estimate the conditional moments of y. If that distribution is known, the integral can be solved using numerical estimation methods such as quadrature. Quadrature is a method of estimating the area under a curve by dividing into many small rectangular intervals and summing the areas of all the intervals (Chapra and Canale, 2006).

Equation 12 computes the probability of a flood with a given magnitude (that is, the probability of a flood of 10,000 cubic feet per second [ft³/s]). Typically, with FFA, the inverse of this function is of interest and the flood magnitude of a given probability is sought (that is, the magnitude of a flood with a 1-percent probability). To estimate the flood magnitude of a specific AEP, optimization can be used to iteratively solve equation 12 across a range of values of y until the AEP of interest is obtained.

When computing the total probability, the conditional probability distribution of the selected climate variable for the year of interest is required. However, assuming the climate variable is the cause of the peak streamflow, nonstationarity means that the climate variable itself is also likely nonstationary, so traditional stationary distribution fitting methods are not applicable to the climate variable of interest. For this study, annual precipitation, temperature, and snow were used as potential explanatory variables for peak streamflow nonstationarity. The conditional distributions for the water year 2020 for climate variables were computed using equations 4–10 with time as the explanatory variable in the regression equation. Lognormal density functions were used for climate variables in equation 12. If future scenarios are of interest, information from large scale climate models could also be used to estimate these climate distributions.

The process used for estimating a flood magnitude for an AEP of interest using a climate-based regression and the law of total probability is outlined below and illustrated in figure 4 using annual precipitation as the climate variable.

Confidence Intervals for Flood Frequency Estimates

Bootstrapped 95-percent confidence intervals were computed to assess the uncertainty in stationary FFA and NFFA estimates for the 1-percent AEP flood. For stationary FFA, the peak streamflow data at each streamgage were sampled with replacement for the same length as the period of record. FFA was then performed on the resampled dataset to obtain a bootstrapped estimate. Several different methods of confidence intervals for nonstationary FFA have been developed (Obeysekera and Salas, 2014). In this study, residual bootstrapping of linear regression equations was used to obtain bootstrapped estimates of the conditional mean, standard deviation, and skews for peak streamflow and climate variable distributions (Davison and Hinkley, 1997; Obeysekera and Salas, 2014). The residuals of the fitted linear regression were sampled with replacement and added to the fitted model estimates to create a new bootstrapped dataset. NFFA was then performed on the bootstrapped dataset. For climate-adjusted NFFA, bootstrapping was performed on the linear regression estimating the conditional climate distribution and the linear regression estimating the conditional peak streamflow. For stationary and nonstationary bootstrap methods, confidence intervals were based on 3,000 bootstrapped estimates and computed as the 2.5-percent and 99.5-percent quantiles of the bootstrapped estimates.

Time-Adjusted Nonstationary Flood Frequency

Time-adjusted NFFA uses time as the only explanatory variable to model peak streamflow nonstationarity. This method is computationally simpler than using a stochastic climate variable and does not require additional data; however, it can be difficult to determine the most appropriate model form.

A time series with a change point can display a variety of statistical characteristics. In some cases, a change point may indicate a step change in the mean or variance of the data. In other cases, there may be a change point within a larger gradual trend, or it may indicate a change in direction of a trend slope. Within a linear regression framework there are several approaches that could be used to model a time series with a change point. Using simple linear regression, a step change in the mean can be modeled using an indicator variable that takes the value of 0 for years equal and prior to the change point and 1 after the change point year. The data before and after the change point can be modeled as a constant or changing in time using an additional explanatory variable for time. A change in slope at a change point can be modeled using an interaction term between a change point indicator and a time variable. Other regression methods such as piecewise regression or local linear regression could also be used in situations where there is a nonlinear pattern of nonstationarity (Chen and others, 2010, Bates and others, 2012), but these methods were not tested in this study.

In cases when statistical tests indicate a trend and a change point, it can be difficult to determine the most appropriate regression model form. For example, figure 5 shows peak streamflow for USGS streamgage 05094000 (South Branch Two Rivers at Lake Bronson, Minnesota), which has a statistically significant change point in water year 1961 (p-value<0.001) and trend (p-value=0.001). Using OLS regression, this peak streamflow time series could be modeled using either a trend or a change point. However, predicted flood magnitudes from these two model specifications were very different in a time-adjusted NFFA. In this case, the change point model estimates a conditional mean of 1,949 ft³/s and a 1-percent AEP streamflow magnitude of 5,287 ft³/s for the year 2020, whereas a trend model estimates a conditional mean of 2,885 ft³/s and a 1-percent AEP streamflow magnitude of 9,462.

Because trend tests can be affected by the presence of change points, it is generally recommended that change points be identified first and that trend tests are performed separately on subsets of the data before and after the change point (Villarini and others, 2009). Site specific exploratory data analysis of hydroclimatic, land use change, and regulation information, as well as regional nonstationarity analyses such as those found in Marti and others (2024), can be used to support the choice of the regression equation form in cases where nonstationarity test results are ambiguous. Examples 1 and 2, discussed in the upcoming sections, demonstrate time-adjusted NFFA at USGS streamgages with a trend and change point in peak streamflows, respectively. For both examples, climate and streamflow data and analyses in Marti and others (2024) were used to characterize the hydroclimate regime and identify probable causal mechanisms for the peak streamflow nonstationarity. Additionally, climate summaries including historical and future scenarios provided by NOAA for each State are used to provide an additional line of evidence to support the use of NFFA (Kunkel and others, 2022).

Example 1

This example estimates the time-adjusted 1-percent AEP flood at USGS streamgage 04108600 (Rabbit River near Hopkins, Michigan) on the western side of the Lower Peninsula of Michigan. This streamgage has 55 peak streamflow observations from water year 1966 through 2020. Mean annual precipitation during the period of record was 37.3 inches and was typically spread throughout the year without a strong seasonal pattern (Marti and others, 2024). Soil is typically near or at saturated conditions during winter and spring (December through May). Peak streamflows during the period of record occur primarily from December through July, with the highest density of peak streamflows in March and April (Marti and others, 2024).

Upward trends in peak streamflow (p-value=0.02, fig. 6A), annual precipitation (p-value<0.01, fig. 6B), and annual temperature (p-value<0.01, not shown) were identified during the period of record. Seasonal trend analyses in Marti and others (2024) showed upward temperature trends in all seasons, with the largest increases in winter (December through February). Upward trends in precipitation were only in winter (December through February) and spring (March through May). Regional analyses using multiple lines of evidence indicate that increases in spring precipitation and changes in winter temperature may be altering snowmelt dynamics and are plausible drivers of changing peak streamflow (Levin, 2024a).

Climate data at this streamgage, as well as regional trend analyses (Levin, 2024a) that show similar precipitation-driven trends in Michigan, provide strong evidence that increases in peak streamflow at this streamgage may be driven by changes in regional climate. Additionally, statewide climate assessments from NOAA predict that increases in precipitation are likely to continue in southern Michigan during the next century, which supports the use of NFFA for this streamgage (Frankson and others, 2022).

Conditional moments of the peak streamflow distribution were computed using a regression equation relating the natural log of peak streamflow (y) and a water year index (wyi), equal to the water year minus 1920. The fitted regression model was y=6.33+0.026×wy_i. Fitted regression coefficients for intercept and slope had p-values of less than 0.001 and 0.01, respectively. The residuals do not show evidence of being heteroscedastic (p-value=0.79) or non-normal (p-value=0.94). The conditional mean for the year 2020 was computed from the regression equation using ω_yi=100. Conditional standard deviation and skew were fit using equations 5 and 6, respectively. The 1-percent AEP quantile was estimated with the method of moments using the quape3 function in the lmomco R package (Asquith, 2023). The conditional and stationary moments and fitted 1-percent AEP quantile are shown in table 1, and figure 7 shows a comparison of stationary and nonstationary estimates across the range of AEPs. The nonstationary estimate of the 1-percent AEP flood for 2020 was 18 percent higher than the stationary estimate.

Table 1.    

Log-Pearson type III distribution moments and 1-percent annual exceedance probability using stationary and time-adjusted flood frequency methods for water year 2020 for U.S. Geological Survey streamgage 04108600 (Rabbit River near Hopkins, Michigan).

Table 1.    Log-Pearson type III distribution moments and 1-percent annual exceedance probability using stationary and time-adjusted flood frequency methods for water year 2020 for U.S. Geological Survey streamgage 04108600 (Rabbit River near Hopkins, Michigan).

Method	Stationary flood frequency	Time-adjusted nonstationary flood frequency
Mean	6.54	6.81
Standard deviation	0.53	0.51
Skew	0.92	0.81
1-percent annual exceedance probability flood estimate	3,345	3,966

Example 2

USGS streamgage 05051600 (Wild Rice River near Rutland, North Dakota) has 58 peak streamflow observations from water years 1960 through 2020. Precipitation has a seasonal pattern with lower precipitation between November and February and higher precipitation in May through July. Median monthly temperatures are typically below freezing from November through March during when streamflow is very low. Streamflow generally increases in the spring as the snowpack melts and precipitation volumes increase. Peak streamflows occur most frequently in the spring primarily owing to snowmelt but can also occur in the summer from large rainstorms (Marti and others, 2024).

Peak streamflows at this location have a statistically significant change point at water year 1992 (p-value<0.001) and a statistically significant upward trend (p-value=0.002, fig. 8A). Mean peak streamflow increased from 82 ft³/s prior to the change point to 468 ft³/s after the change point, which is a 5.7-fold increase. Annual precipitation at this location also has an upward change point in 1990 (p-value=0.016, fig. 8B) and an upward trend (p-value=0.017). An analysis of the timing of peak streamflows and daily streamflows is shown in figure 9A (Marti and others, 2024). In the mid-1990s, there was an abrupt change in the seasonality of peak streamflows. Although peaks prior to the mid-1990s occurred primarily in the spring, after the mid-1990s there was an increase in peaks in May through August (fig. 9A). Additionally, daily streamflows and lower flows increased abruptly in the mid-1990s (fig. 9B). Modelled annual soil storage data indicate that annual soil storage increased concurrent with the change point in precipitation (fig. 9C). The abrupt increase in precipitation and antecedent soil storage during the summer is a likely driver of the change in seasonality and magnitude of peak streamflows.

Climate data for this streamgage, as well as regional trend analyses (Ryberg and others, 2024), provide strong evidence that increases in peak streamflow at this streamgage may be caused by changes in regional climate. Additionally, the concurrent change point in precipitation and soil moisture, as well as the apparent abrupt change in peak streamflow seasonality and daily streamflow at the time of the change point, provide strong evidence that peak streamflows in this location have undergone an abrupt shift rather than a gradual trend and that time-adjusted NFFA analyses may be more accurately represented with a change point model.

Peak streamflows were transformed using the natural log (y) and a change point indicator variable (cpind) was computed that was equal to 0 for years 1960 through 1992 and equal to 1 from 1993 to 2020. The period after the change point does not have a statistically significant trend (p-value=0.80) and the regression equation used only the change point indicator as an independent variable. The fitted regression equation was y=4.40+1.74×cpind. Fitted regression coefficients for intercept and change point indicator had p-values of less than 0.001. Residuals from this regression did not show significant heteroscedasticity (p-value=0.062) or non-normality (p-value=0.104). The conditional mean after the change point is estimated by the regression equation as ${\mu }_{y|cpind}$=4.4+1.74×1=6.14. Conditional standard deviation and skew were derived from the fitted regression residuals and standardized time series using equations 6 and 9, respectively. The 1-percent AEP quantile was estimated with the method of moments using the quape3 function in the lmomco R package (Asquith, 2023). The conditional and stationary moments and fitted 1-percent AEP quantile are shown in table 2, and figure 10 shows a comparison of stationary and time-adjusted nonstationary estimates across the range of AEPs.

Table 2.    

Log-Pearson type III distribution moments and 1-percent annual exceedance probability using stationary and nonstationary methods for water year 2020 flood frequency methods for U.S. Geological Survey streamgage 05051600 (Wild Rice River near Rutland, North Dakota).

Table 2.    Log-Pearson type III distribution moments and 1-percent annual exceedance probability using stationary and nonstationary methods for water year 2020 flood frequency methods for U.S. Geological Survey streamgage 05051600 (Wild Rice River near Rutland, North Dakota).

Method	Stationary flood frequency	Time-adjusted nonstationary flood frequency
Mean	5.25	6.15
Standard deviation	1.65	1.4
Skew	−1.00	−0.8
1-percent annual exceedance probability flood estimate	2,586	5,321

Climate-Adjusted Nonstationary Flood Frequency Analysis

Climate-adjusted NFFA uses a climate variable to model the nonstationarity in peak streamflow. This method requires a regression equation relating peak streamflow to the climate variable and the probability distribution of the climate variable for the year of interest. Climate-adjusted NFFA is more computationally complex than time-adjusted NFFA and there may be additional uncertainty from the inclusion of the estimated climate probability distribution.

Regional trend attribution studies have identified large-scale changes in climate as a driver of peak streamflow nonstationarity in the study area (Levin and Holtschlag, 2022; Sando and others, 2022). Changes in precipitation magnitude and timing have been identified as likely causes of increasing peak flows in the upper Midwest (Levin and Holtschlag, 2022; Sando and others, 2022). Decreasing peak flows have been attributed to land use change, groundwater withdrawals, decreases in soil moisture, increased temperature, and decreases in snowpack (Neri and others, 2019; Sando and others, 2022). For this study, only annual precipitation, temperature, and snowfall were considered as explanatory variables for regressions estimating peak streamflow. These variables were chosen because of their wide-scale data availability across the region and because they have been used in previous related publications characterizing nonstationarity in this region (Ryberg, 2024).

Because the conditional variance and skew for NFFA analyses are computed from regression residuals, it is important that fitted regressions have residuals that are homoscedastic and normally distributed. Failure to meet these assumptions may result in biased estimates of these distribution parameters. The assumptions of homoscedasticity and normality were tested with the Breusch-Pagan test for homoscedasticity (Breusch and Pagan, 1979) and the probability plot correlation test (Vogel, 1986), respectively. Time-adjusted conditional distribution moments for climate variables were computed for the year 2020 using equations 5–11. Annual climate variables were fit to lognormal probability distributions. The conditional flood and the conditional climate variable distribution were used in the total probability equation (eq. 12) and optimization was used to determine the streamflow magnitude for the 1-percent AEP flood magnitude.

At some streamgages with change points in peak streamflows, there was a change in the relation between peak streamflow and the explanatory climate variable before and after the change point. This change can happen when a change in temperature or soil moisture modifies the relation between precipitation and streamflow (Woodhouse and others, 2016). For example, a step increase in precipitation can increase soil moisture throughout the year. The change in antecedent soil moisture can cause greater streamflow for a given precipitation event after the change point year. In these cases, a change point indicator variable equal to 0 for observations prior to and including the year of the change point and 1 after the change point can be included in the regression equation between peak streamflow and the climate variable to improve the fitted regression equation (refer to ”Example 4” section).

The following examples demonstrate the process of climate-adjusted NFFA at two streamgages. The examples estimate the 1-percent AEP flood using climate conditions for the year 2020.

Example 3

USGS streamgage 05107500 (Roseau River at Ross, Minn.) has 80 peak streamflow observations from water year 1929 through 2020. This streamgage is located along the northern border of Minnesota and streamflow is strongly influenced by snowmelt. Precipitation has a seasonal pattern with lowest precipitation from November through March and highest precipitation in May through July. Median monthly temperatures are typically below freezing from November through March during which time streamflow is very low. Streamflows increase in the spring as the snowpack melts and precipitation volumes increase (Marti and others, 2024).

Peak streamflows at this location have a marginally significant upward change point at water year 1961 (p-value=0.08) and a statistically significant upward trend (p-value=0.02, fig. 11A). Annual precipitation at this location has an upward change point (p-value=0.007) and an upward trend (p-value=0.002, fig. 11B) driven largely by upward trends in winter and spring (Marti and others, 2024). Temperature at this location also increased and had a statistically significant change point (p-value<0.001) in 1979. Increased precipitation accompanied by warmer temperatures, particularly during the winter, which can cause rain-on-snow events and changes in snowmelt dynamics, are consistent as a causal mechanism for the increase in peak streamflows at this location. Climate data for this streamgage, as well as regional trend analyses (Marti and others, 2024), provide strong evidence that increases in peak streamflow at this streamgage may be caused by changes in regional climate. Statewide climate assessments by NOAA predict that high precipitation rates will continue in this region through this century, supporting the use of NFFA at this streamgage (Runkle and others, 2022).

Conditional moments for peak streamflows were based on the relation between peak streamflows and annual precipitation (fig. 12). The fitted regression equation between the natural log of peak streamflows (y) and natural log of annual precipitation (w) was y=1.78+1.84×w. Residuals for this regression equation do not show evidence of heteroscedasticity (p-value=0.41) or non-normality (p-value=0.21). Conditional moments for peak streamflow were computed using equations 5–8.

Conditional moments of annual precipitation for the year 2020 were computed from the regression equation relating the natural logarithm of annual precipitation (w) and a change point indicator (cpind) equal to zero prior to and including water year 1961 and equal to 1 after that: w=2.98+0.14×cpind. Residuals from this equation do not show evidence of heteroscedasticity (p-value=0.41) or non-normality (p-value=0.21). The conditional mean of log-transformed annual precipitation for the year 2020 was 3.12, and the standard deviation was 0.179. The conditional probability density function for the log-transformed annual precipitation for the year 2020 was estimated with the dnorm function in R.

Equation 12 was computed across the range of precipitation values from the 0.0001 to 0.9999 quantiles of the computed conditional precipitation distribution. Optimization of equation 12 was performed in R with the optim() function, starting with an initial flood estimate equal to the minimum observed peak streamflow and iterating until the estimated annual exceedance probability was within 0.0002 of the target AEP value (in this case, 1 percent). Estimates of the 1-percent AEP flood using stationary, time-adjusted NFFA, and climate-adjusted NFFA methods are shown in table 3. The conditional mean is not reported in table 3 for climate-adjusted NFFA because there are a range of conditional means used in equation 12. Plots of the annual exceedance probability across a range of values using stationary, time-adjusted NFFA, and climate-adjusted NFFA methods are shown in figure 13.

Table 3.    

Log-Pearson type III distribution moments and 1-percent annual exceedance probability using stationary, time-adjusted nonstationary, and climate-adjusted nonstationary flood frequency methods for water year 2020 for U.S. Geological Survey streamgage 05107500 (Roseau River at Ross, Minnesota).

[--, no data]

Table 3.    Log-Pearson type III distribution moments and 1-percent annual exceedance probability using stationary, time-adjusted nonstationary, and climate-adjusted nonstationary flood frequency methods for water year 2020 for U.S. Geological Survey streamgage 05107500 (Roseau River at Ross, Minnesota).

Method	Stationary	Time-adjusted nonstationary flood frequency	Climate-adjusted nonstationary flood frequency
Mean	7.43	7.73	--
Standard deviation	0.7	0.68	0.61
Skew	−0.55	−0.55	−0.38
1-percent annual exceedance probability flood estimate	6,486	8,463	7,993

Example 4

USGS streamgage 06326500 (Powder River near Locate, Montana) has 78 peak streamflow observations from water year 1938 through 2020. Precipitation has a seasonal pattern with lowest precipitation from November and March and highest precipitation in May and June. Median monthly temperatures are typically below freezing from November through March during which time streamflow is very low. Similar to other northern locations in the study area, peak streamflow at this location is influenced by snowfall and snowmelt, which can persist until May or June (Marti and others, 2024).

Peak streamflow at this location has a downward change point in 1972 (p-value=0.001, fig. 14A). There was no trend or change point in precipitation at this location, but nonstationarity was detected in annual snow and annual mean temperature (fig. 14). Annual mean temperature had an upward change point in 1986 (p-value<0.001) and annual snowfall had downward change point in 1979 (p-value=0.02). Previous studies have attributed downward trends in peak streamflow in Montana to combination of precipitation and temperature effects that are difficult to separate (Sando and others, 2022). Nonstationarity in temperature can modify the relation between streamflow and precipitation or snowfall by changing evapotranspiration, soil moisture, and snowmelt processes (Woodhouse and others, 2016). A decrease in snowfall is a likely driver of downward changes in peak streamflow; however, an apparent change in the relation between peak streamflow and snowfall before and after the change point is likely due to the change in temperature (fig. 15).

Conditional moments for peak streamflow were based on the relation among log transformed peak streamflow (y), annual snowfall (w), and a change point indicator variable (cpind) equal to 0 prior to and including the change point year and equal to 1 after (fig. 15). The fitted regression equation was y=7.35+1.74×w–0.56×cpind. All fitted coefficients were statistically significant with p-values<0.01. Residuals for this equation did not show significant heteroscedasticity (p-value=0.54) or non-normality (p-value=0.18). Conditional moments were computed using equations 5–8 (table 4).

Conditional moments for annual snowfall were computed for the year 2020 using the relation between the natural logarithm of annual snowfall (w) and time, using a change point indicator (cpind) equal to zero prior to or equal to 1979 and 1 afterwards, with the equation w=1.11–0.16×cpind. The conditional mean of log-transformed annual snowfall for the year 2020 of 0.95, and the standard deviation was 0.228. The conditional probability density function for the log-transformed annual snowfall for the year 2020 was estimated with the dnorm function in R. Equation 12 was computed across the range of snowfall distribution. Optimization of equation 12 was performed in R with the optim() function, starting with an initial flood estimate equal to the minimum observed peak streamflow and iterating until the estimated annual exceedance probability was within 0.0002 of the target AEP value (in this case, 1 percent).

Estimates of the 1-percent AEP flood using stationary, time-adjusted NFFA, and snowfall-adjusted NFFA methods are shown in table 4. The conditional mean is not reported for climate-adjusted NFFA because there is a range of conditional means used in equation 12. Plots of the annual exceedance probability across a range of values using stationary, time-adjusted NFFA, and climate-adjusted NFFA methods are shown in figure 16. In this case, time-adjusted and climate-adjusted methods produce similar results, both of which estimate lower flood magnitudes than the stationary method.

Table 4.    

Log-Pearson type III distribution moments and 1-percent annual exceedance probability using stationary, time-adjusted nonstationary, and climate-adjusted nonstationary flood frequency methods for water year 2020 for U.S. Geological Survey streamgage 06326500 (Powder River near Locate, Montana).

[--, no data]

Table 4.    Log-Pearson type III distribution moments and 1-percent annual exceedance probability using stationary, time-adjusted nonstationary, and climate-adjusted nonstationary flood frequency methods for water year 2020 for U.S. Geological Survey streamgage 06326500 (Powder River near Locate, Montana).

Method	Stationary	Time-adjusted nonstationary flood frequency	Climate-adjusted nonstationary flood frequency
Mean	8.81	8.52	--
Standard deviation	0.84	0.77	0.56
Skew	−0.28	−0.43	−0.32
1-percent annual exceedance probability flood estimate	39,608	23,310	24,953

Regional Time-Adjusted Nonstationary Flood Frequency Analyses

NFFA using time as the explanatory variable was performed at candidate streamgages in the study area. There were 107 streamgages that had statistically significant change points as well as significant trends, similar to example 2. It was beyond the scope of this study to do detailed site-specific data analyses at all of these streamgages to determine which model form would be more appropriate. Instead, if a change point and trend were identified at a location, it was modeled as a change point as is suggested by Villarini and others (2009).

The model form (change point or trend) of the regression equation used for time-adjusted NFFA and how well the regression residuals met the assumptions of OLS regression are shown in figure 17. Of the 153 candidate streamgages, 103 were modeled with a change point indicator model and 34 were modeled as trends. There were 16 streamgages for which time-adjusted NFFA was not performed because fitted regression residuals exhibited extensive heteroscedasticity or non-normality. Of the 137 streamgages for which time-adjusted NFFA was performed, regression models at 101 streamgages had residuals that did not show significant heteroscedasticity or non-normality or whose residuals could be adequately modeled using a secondary error model (eq. 10), making them good candidates for this NFFA method (fig. 17). The other 36 streamgages had regression equations for which the regression assumptions were marginally met.

The percent difference in the stationary and time-adjusted nonstationary 1-percent AEP flood estimate was examined at candidate streamgages in the study area (fig. 18). Comparisons of FFA and NFFA methods in this report should be taken only as a screening-level analysis to give an indication of the general magnitude of change. Because it was outside the scope of this study to use the expected moments algorithm (EMA) method of FFA as recommended by Bulletin 17C, the stationary FFA methods used in this study may not match currently published estimates. Despite these limitations, these comparisons give a general estimate of the change in flood magnitudes owing to nonstationarity and may identify regional patterns.

Differences in flood estimates from time-adjusted NFFA and stationary FFA across the region ranged from −80 percent to 100 percent. The greatest percent differences between stationary FFA and time-adjusted NFFA estimates for the 1-percent AEP flood magnitude were in the northern part of the study area. Upward changes were greatest in eastern North Dakota and Minnesota, whereas the greatest downward changes in estimates were in western North Dakota, South Dakota, and eastern Montana (fig. 18). Differences in the 1-percent AEP were less than 20 percent at most streamgages in Missouri, Michigan, and Illinois.

Time-adjusted NFFA using linear regression was applicable throughout most of the study area. The assumptions of linear regression were adequately met for trend and change point models at most streamgages. Although most streamgages were fit with change point models, nearly all the streamgages with change points also had statistically significant trends. In some cases, site-specific exploratory analysis at individual streamgages (examples 1–4) may lend support for different modeling forms than those selected in this report.

Regional Climate-Adjusted Nonstationary Flood Frequency

Climate-adjusted NFFA was performed at streamgages with a nonstationary annual climate variable that exhibited a statistically significant trend or change point at the 95-percent confidence level that was consistent with the direction of significant change in the peak streamflows (fig. 2). The climate variable selected for use in NFFA at streamgages across the study area is shown in figure 19. Regression equations for upward trends and change points in peak streamflow at streamgages throughout the central part of the study area were fit using annual precipitation if a concurrent upward trend or change point was identified. Regression equations for downward peak streamflow trends and change points in the western part of the study area and in Wisconsin and upper Michigan were fit using either annual temperature or annual snowfall. If temperature and snowfall were nonstationary at a streamgage, the variable with the greater regression coefficient of determination (R²) was chosen. The R² measures the proportion of the variability in peak streamflow that is explained by the explanatory variable (annual temperature or snowfall). Of the 153 candidate streamgages, 70 were predicted using annual precipitation, 21 were predicted using annual mean temperature, and 7 were predicted using annual snowfall. There were 55 streamgages for which a climate variable could not be identified as a causal mechanism for peak streamflow under the criteria used in this project.

Of the 98 streamgages for which an appropriate climate variable was identified, 65 had fitted regression models with residuals that did not show significant heteroscedasticity or non-normality or whose residuals could be adequately modeled using a secondary error model (eq. 10), making them good candidates for this NFFA method. An additional 33 streamgages had regression equations for which the regression assumptions were marginally met. Climate-adjusted NFFA analyses at these streamgages can benefit from closer examination and may need different explanatory climate variables, a more robust model form, or an alternative NFFA method.

The percent difference between stationary and climate-adjusted NFFA estimates of the 1-percent AEP flood is shown in figure 20. Results were similar to those from the time-adjusted analyses, with the greatest upward changes in eastern North and South Dakota and Minnesota. A comparison of climate and time-adjusted NFFA estimates at streamgages where both analyses were performed is shown in figure 21. Differences between climate- and time-adjusted estimates agreed within 25 percent of each other in about 75 percent of the cases.

Climate analyses were applicable in about 60 percent of the streamgages with observed peak streamflow nonstationarity. Many streamgages were excluded from analysis because none of the three annual climate variables could be identified as a probable cause of nonstationarity, either because they were stationary or because the direction of the climate trend or change point was not consistent with a causal mechanism for the peak streamflow trend or change point. In cases where a climate variable was identified as a likely driver of nonstationarity, the method produced 1-percent AEP estimates that were similar in magnitude to those estimated by time-adjusted NFFA methods. Additionally, regression equations using climate variables generally explained a larger portion of the variance than time-based regression equations, with mean coefficient of determination of 0.34 for climate-based regressions and 0.13 for time-based regression equations.

Uncertainty in FFA is high, even under stationary conditions. Sources of uncertainty in stationary FFA include sampling error, uncertainty in estimated sample moments, and the fit of the selected distribution. Many of the same sources of uncertainty in stationary FFA also apply to NFFA. Additionally, NFFA methods in this study are affected by uncertainty from fitted regression equations. Although a full examination of the controls on uncertainty in regression-based NFFA are beyond the scope of this study, the same factors that affect uncertainty in regression are likely to affect NFFA estimates, including sample size, error variance of the fitted regression, and the presence of outliers or influential data points.

Median lower and upper bootstrapped 95-percent confidence intervals were −28 and +38 percent of the estimated value for stationary FFA, −32 and +43 percent for time-adjusted NFFA, and −26 and +55 percent for climate-adjusted NFFA (fig. 22). At some streamgages, NFFA confidence intervals were excessively large, with upper confidence intervals greater than 400 percent. Large confidence intervals occurred at streamgages with shorter periods of record (less than 70 years), those that had large regression error variance, or several large outliers in the data. In general, bootstrapped confidence intervals were larger at streamgages with shorter periods of record where regression slopes were more sensitive to variability in the bootstrapped data and produced a wider range of conditional means. The coefficient of variability of bootstrapped conditional means was twice as large in streamgages with less than 70 years of record than it was at sites with longer periods of record (coefficient of variability of 2.5 and 1.1, respectively). The large range of conditional means reflects the sensitivity of OLS regression slopes to data variability and outliers, particularly when samples sizes are smaller. Additional research can help characterize the factors affecting uncertainty in NFFA estimates; however, this preliminary examination indicates that regression equations that explain a greater portion of the variability and larger sample sizes may produce flood estimates with less uncertainty.

Limitations

Analyses in this study should be regarded as preliminary and do not supersede other published flood frequency estimates. There are several limitations of analyses performed in this study. Only annual precipitation, snowfall, and temperature were considered as potential explanatory variables in this study. These explanatory variables were chosen because of their availability across the large study area and their previous use in characterization of nonstationarity in the study area (Ryberg, 2024); however, subannual climate variables may be better predictors of flood magnitudes at some locations. Additionally, these three climate variables do not represent the full range of flood-generating mechanisms in the region. Changes in land use from agriculture, urbanization, and tile drainage may also have effects on flood magnitude but were beyond the scope of this study. Additionally, only linear OLS regression was considered for estimating the relation between explanatory variables and peak streamflow. Other regression methods such as piecewise regression, robust regression, or multiple regression, which considers the mutual effects of more than one explanatory variable were not considered but may result in better model fits and reduce the uncertainty in estimates. Finally, site-specific exploratory data analyses (examples 1–4) were not completed for the regional application of the method. These regional analyses were performed as a preliminary assessment of the applicability of the method and to identify areas where changes were greatest in the estimated 1-percent AEP. Closer inspection of climate and streamflow trends at individual sites may support the use of different explanatory variables or model forms than those used in this study.