Regression modeling of combined sewer overflows to assess system performance

Climate changes including more frequent and extreme wet weather events (Easterling et al. 2017) stress aging wastewater and stormwater infrastructure systems designed for the climate of the past (Rudberg et al. 2012; Schoen et al. 2015). These stresses are particularly acute in combined sewer systems (CSS), in which stormwater is collected in the same pipes as domestic and/or industrial wastewater. Thousands of CSOs occur each year in the United States (US Environmental Protection Agency 2022a), Europe (Environmental Agency 2020; Scottish Water 2020), and elsewhere around the world. Changing precipitation patterns may lead to more severe and more frequent incidents of combined sewer overflow (CSO), releases of untreated raw sewage mixed with stormwater from CSS (Astaraie-Imani et al. 2012; Fortier & Mailhot 2015; Mahaut & Andrieu 2019; Botturi et al. 2021). Though climate change is forecast to reduce rainfall events and wet extremes in some areas (Christidis & Stott 2022), much of North America (Easterling et al. 2017) and northern Europe (Christidis & Stott 2022) are expected to face more rainfall and worsening wet extremes. In these regions, urbanizing areas with combined sewers are particularly vulnerable, as increasing urbanization, absent mitigation of additional flow volumes, also increases the frequency and volume of CSO discharges (Semadeni-Davies et al. 2008; Abdellatif et al. 2015). Indeed, emerging research on predicting CSO occurrence suggests CSO events will worsen in urban areas under climate change (Semadeni-Davies et al. 2008; Schroeder et al. 2011; Astaraie-Imani et al. 2012; Sandoval et al. 2013; Fortier & Mailhot 2015; Yu et al. 2018).

When CSOs harm a community or surrounding environment, government regulators may respond by taking enforcement actions requiring CSS municipalities to reduce CSO discharges or otherwise act to improve water quality. In the United States, this typically takes the form of a consent decree from the US Environmental Protection Agency (USEPA) ordering the reduction or elimination of CSO discharge volumes (US Environmental Protection Agency 1994). Some U.S. municipalities end CSO discharges altogether by separating their sewer systems from their stormwater collection systems, while other municipalities, particularly those who serve historically marginalized communities, cannot afford to separate their systems. Existing methods for assessing CSS performance relative to planned CSO reductions rely upon complex hydraulic modeling approaches that are not easily transferable and yield results with wide uncertainty ranges (Andrés-Doménech et al. 2012; Bachmann-Machnik et al. 2021; Roseboro et al. 2021). Yet, simpler approaches to CSS assessment with lower uncertainties often have trade-offs, such as ignoring CSO magnitudes entirely and instead measuring only water quality (Madoux-Humery et al. 2013; Romei et al. 2021). To understand the effects of climate change on CSO discharges, rather than measuring water quality, it is simpler to assess the frequency and volume of CSO discharges. Using the city of Cumberland, Maryland as a case study, we build on prior research to develop a simple yet transferable statistical model to predict the occurrence and volume of CSO incidents and propose a new performance indicator of critical intensity (I_crit) based on CSO occurrence and rainfall intensity to assess CSS performance.

Schroeder et al. (2011) proposed a simple statistical approach to predict whether a CSO will occur from precipitation data. They identified a critical rainfall depth (N_crit), the total rainfall amount (in mm) of an antecedent rainfall event, above which a CSO was probable (≥85%) to occur and below which a CSO was improbable (≤17%) to occur, for each of four watersheds in Berlin, Germany, in 2007. Expanding on this work, Fortier & Mailhot (2015) identified N_crit and fit a sigmoid probability function for CSO incidence on the antecedent rainfall depth for outfalls in Quebec, Canada. Both Schroeder et al. (2011) and Fortier & Mailhot (2015) used antecedent rainfall depth (mm of precipitation) in their modeling procedures. Consistent with these, Sandoval et al. (2013) found that rainfall depth is an effective predictor of CSO volume and duration. Sandoval et al. (2013) also found that the maximum precipitation intensity of the antecedent rainfall event is the most effective predictor of CSO volume, duration, and flow rate. More recently, Yu et al. (2018) and Hung et al. (2020) suggested that the intensity (mm/h), whether recorded at intervals of 1 min (Yu et al. 2018) or calculated as an average (Hung et al. 2020), of an antecedent rainfall event is more effective at identifying probabilistic CSO thresholds in outfalls in urban watersheds, and Li et al. (2021) demonstrated the use of both depth and intensity of rainfall to effectively predict the volumes and peak rates of flows within CSS pipes. Underscoring all of this, Todeschini (2012) found that, as the climate changes, increasing rainfall intensities present the most significant challenges for urban stormwater management, as more intense average rainfall events will increase CSO discharge volumes regardless of other changes in annual precipitation patterns.

Here we demonstrate a simple and transferable statistical approach using precipitation to predict CSOs, which can be used to assess CSS performance over time. The case study for this research is Cumberland, Maryland, which operates a small wastewater treatment plant (average flow 8–13 million gallons (∼36.4–49.2 million L) per day with a treatment capacity of 15 million gallons (∼56.8 million L) per day) (US Environmental Protection Agency 2022b) served by a small CSS that includes a moderately environmentally disadvantaged, low-income urban area, and that sits mostly within the Wills Creek and Potomac Lower North watersheds. The USEPA requires Cumberland to reduce its annual CSO volumes to specific targets by 2023. While other Maryland CSS communities separated their stormwater and sewer systems, Cumberland has adopted an approach to reduce CSOs by rehabilitating the existing system, implementing asset management, and expanding the carrying and storage capacity of the CSS. Cumberland has documented all CSO incidents since 2005 and made that information publicly available, in accordance with 2005 Maryland state regulations on tracking sewage spills (Maryland Department of the Environment 2005). To develop our predictive model of CSO incidence and volume from precipitation variables, we analyze 16 years of Cumberland CSO reports, which include data on overflow incidence, overflow causes, discharged volumes, and times and locations of overflows, spanning January 2005 through December 2020, in both the Wills Creek and Potomac Lower North watersheds, in combination with historical precipitation data from Phase 2 of the North American Land Data Assimilation System (NLDAS-2) covering the same period. We assume no change in urbanization and impervious cover over the time period because urbanization and population are not expanding in Cumberland.

North American Land Data Assimilation System Phase 2 (NLDAS-2) precipitation data at 1/8th-degree geospatial resolution and hourly temporal resolution were used for analysis. NLDAS hourly precipitation data are developed using atmospheric forcing and routing models from the execution of four Land Surface Models (LSMs) on a data set assimilated from the rainfall data of the North American Regional Reanalysis (NARR) and Climate Prediction Center (CPC), both of which are based on direct rainfall gauge observations from across the continent that have been disaggregated and adjusted for topography using the Parameter-elevation Regression on Independent Slopes Model (PRISM) (Xia et al. 2009). Due to the coarse geospatial resolution, only one set of hourly precipitation data was available for the Cumberland area. Discrete rainfall events were identified as rainfall periods separated by either 12 + consecutive hours of precipitation intensity below 0.2 mm/h (Fortier & Mailhot 2015), or 2 + consecutive hours of zero precipitation (Hung et al. 2020). The latter criterion was applied when 12 + consecutive hours of precipitation intensity below 0.2 mm/h failed to discretize rainfall events enough to pair individual rainfall events with daily CSO events.

We adopted the approach of Fortier & Mailhot (2015) to pair CSOs with rainfall events. Each CSO discharge was assigned a ‘start time’ equal to the discovery time in the CSO incident report, if such a discovery time was given, or 12.00 noon when no discovery time was given. Forty eight percent of CSOs included a discovery time, leaving 52% with a 12.00 noon start time assigned. For each CSO incident, we selected discrete rainfall events for pairing if they occurred between 48 hours before the ‘start time’ and 12 hours after the start time (considering the latter permissible to accommodate both uncertainty in CSO reports and the coarse geospatial resolution of the NLDAS-2 rainfall data). When a CSO incident had multiple discrete rainfall events that could be paired to it, one event was chosen based on greater total rainfall depth, greater rainfall intensity, and occurrence of the day before the CSO date, in that order of priority. This procedure and prioritization are adapted from the methodology of pairing CSO and rainfall events used by Fortier & Mailhot (2015). Additionally, rainfall amounts from any other discrete rainfall events between this selected event and the CSO ‘start time’ were then added into the selected event to create a composite total picture of antecedent precipitation for each CSO incident. This process resulted in the tabulation of 5,405 rainfall events with and without CSO incidence, including 2,696 in the Wills Creek watershed and 2,709 in the Potomac Lower North watershed, with each rainfall event's maximum and average intensities, antecedent dry period, and total rainfall depth recorded from the NLDAS-2 data set. The associated CSOs were grouped by watershed and date, and then each group was summed into a single CSO event, to model CSOs as watershed-scale responses to wet weather. This resulted in consolidated lists of 715 paired CSO-rainfall events in the Wills Creek watershed and 599 paired CSO-rainfall events in the Potomac Lower North watershed, which were used for all subsequent analysis.

Using the CSO-rainfall event data, we modeled CSO incidence as a response variable to precipitation conditions (predictor variables) to identify and evaluate CSO incidence thresholds for the Potomac Lower North and Wills Creek watersheds. First, we identified a critical total rainfall depth N_crit (mm) which best discriminates CSO events from non-CSO events for each year of data in each watershed, following the method demonstrated by Schroeder et al. (2011), i.e. maximizing the average of the proportions of CSO-rainfall events above the threshold and of non-CSO rainfall events below the threshold. We performed logistic regressions using R (R Core Team 2021) of CSO incidence on each of rainfall depth (P, mm), maximum rainfall intensity (I_max, mm/h), average rainfall intensity (I_avg, mm/h), and antecedent consecutive dry period (CD, h), to identify the best performing predictors. Event maximum intensity was identified as an effective predictor of CSO incidence based on logistic regressions producing higher odds ratios without vast uncertainties, and as such, we identified a critical maximum intensity I_crit (mm/h) that best discriminates CSO events from non-CSO events for each watershed in each year, adapting the method of Schroeder et al. (2011) to use event maximum intensity instead of depth. The yearly sets of I_crit and N_crit values were comparatively evaluated for their ability to predict CSOs, and for correlations to year-to-year climate variables such as total annual rainfall, counts of all rainfall events, and counts of more severe rainfall events, derived from the NLDAS-2 data.

Using the CSO-rainfall event data, we performed ordinary least-squares (OLS) regression for gross CSO volume on associated event rainfall depth (P) for the CSO events in each watershed and year, based on Fortier & Mailhot (2015). Further analysis, performed in R, of bivariate correlations and regression diagnostics found that while the event rainfall depth is correlated with CSO volume, OLS regression produced highly skewed residuals (see Supplementary Material, Figure S1). To correct this, we transformed the data using natural logarithms, and then performed regressions on the transformed (lognormal) data. We used R programming to develop, evaluate, and compare multiple regression models of CSO volume (lnV) on the natural logarithms of candidate predictors including event rainfall depth (lnP), maximum rainfall intensity (lnI_max), average rainfall intensity (lnI_avg), and interactions between these predictors, using the full 16-year data sets for all CSO incidents in each of the Wills Creek and Potomac Lower North watersheds. We then applied the two most promising candidate models to the single-year CSO-rainfall event data sets to judge their historic consistency in correlation and diagnostics to identify the best reliable regression model. This model can be applied over each year (or other designated periods of time) to evaluate CSS performance.

CSO incidence is treated as a binary response variable to precipitation conditions, with total rainfall depth, maximum intensity, average intensity, and length of the dry period preceding rainfall as possible predictor variables. Welch two-sample t-tests were performed on the means of each variable in rainfall events with CSOs versus those without CSOs, and found statistically significant differences in these means (p < <0.001) in each watershed (see Supplementary Material, Figure S2).

N_crit values ranged between 3.1 and 6 mm in the Wills Creek watershed, and between 3.3 mm and 7.0 mm in the Potomac Lower North watershed, for each year of data from 2005 through 2020. In any given year, N_crit values correctly discriminate between 72 and 87% of rainfall events. These probabilities were lower than expected and generally lower than the values demonstrated in Schroeder et al. (2011).

Importantly, in year-to-year summary data with the annual sum of NLDAS-2 rainfall events representing annual rainfall totals, we found that N_crit values in each watershed positively correlate to the total annual rainfall recorded over Cumberland (Pearson's r = 0.626 (p = 0.0095) in the Wills Creek watershed; Pearson's r = 0.535 (p = 0.0329) in the Potomac Lower North watershed). Because N_crit is not independent of yearly precipitation changes, N_crit is not a good indicator of system performance over time.

On the yearly scale, odds ratios for logistic regressions of CSO incidence on average intensity ranged from 1.38 to 8.47, as seen in Figure 3. This was consistently the highest odds ratio in any given year, but with great variability. Also, at the annual scale, much higher odds ratios also had much higher uncertainties in the odds ratio. Having odds ratios with high uncertainties and variable magnitudes meant that average intensity would not be useful for assessing system performance. Yearly logistic regressions of CSO incidence on maximum intensity produced more consistent yearly odds ratios (ranging from 1.37 to 2.91), with moderate uncertainties. Odds ratios on maximum intensity were also statistically significantly higher than those for CSO incidence versus total depth in each watershed in each year. Therefore, we used maximum intensity to define a new probabilistic CSO threshold, the critical intensity I_crit, in each watershed and year.

CSO volume is treated as a response to the related rainfall event's precipitation depth (P, mm), average intensity (I_avg, mm/h), maximum intensity (I_max, mm/h), and antecedent consecutive dry period (CD, h). These variables have generally right-skewed distributions with high variance and spread, and moderate, manageable amounts of correlation between them. Correlation matrices run over the full data set in both watersheds found P to have the strongest correlation with CSO volume (V, gal, where 1 gal = 3.7854 L) at Pearson's r = 0.24 in the Wills Creek watershed and Pearson's r = 0.38 in the Potomac Lower North watershed (Supplementary Material, Table S1). However, OLS regression of V on P produced invalid regression results; we observed that the residuals of these regressions were highly skewed, and there were many outliers (see Supplementary Material, Figure S1). These issues were resolved by treating all variables as lognormal distributions. In turn, the univariate distributions of lnV, lnP, lnI_avg, and lnI_max are approximately normal. Bivariate linear regressions explained highly inconsistent amounts of variance in lnV that were not always statistically significant against a null hypothesis of no correlation (with α = 0.05), and often with fit that still indicated nonlinearity. Thus, we deemed multiple linear regression to be a preferred approach.

This multiple regression model performs better than the simpler regressions of the natural logarithm of volume (lnV) on the natural logarithm of either event precipitation depth (lnP) or intensity (lnI_avg) alone. Figure 6 presents how the yearly coefficients of determination compare to the bivariate regression of lnV vs lnP, showing the bivariate correlation to be much less reliable. Although the bivariate correlation is sometimes higher than that of the multiple regression, the multiple regression's intensity terms prevent such drastic losses of correlation as can be sometimes observed in the bivariate model, and dramatically reduce the variance among R² values of the model. Thus, the multiple regression is the preferred model.

The consistency of this model from year to year makes it a reasonable choice to assess or identify changes in CSS performance, including to evaluate the impacts of climate change and interventions intended to build resilience to climate change. For example, in 2021, Cumberland installed a 50 million gallon (∼227 million L) capacity storage volume to mitigate CSOs. The city plans to further increase capacity and mitigate CSOs by installing 78 in (1.98 m)-diameter overflow carrier pipes in the mid-2020s. Changes in both the I_crit thresholds and the regression results measured over each year (or a longer time period) could be used to evaluate the success of these interventions in mitigating CSOs in Cumberland, or adapted for analysis in other communities with comparable data availability.

Using data on combined sewer overflows and hourly precipitation in Cumberland, Maryland over the years 2005–2020, we identified annual thresholds of 1-hour precipitation intensity above which CSO incidence is predicted (and below which it is not predicted) in each year. Then, we built a regression model to predict CSO volume from the precipitation depth and average intensity of the preceding rainfall event. These annual thresholds and the regression model can be used in tandem to assess the performance of a combined sewer system (CSS), such as Cumberland's, over time.

We developed this as a two-pronged model to predict CSO incidence and volume. Changes to the model outputs can be read over individual time periods (of 1 year or longer) as indicators of changing CSS performance over time. The model consists of:

• 1.

  Identifying a critical maximum rainfall intensity (I_crit) which most successfully discriminates CSO-associated rainfall events from non-CSO rainfall events over time

• 2.

  Performing multiple regression of the natural logarithm of CSO volume (gal) on the natural logarithms of the associated rainfall event's total precipitation depth (mm) and average intensity (mm/h) with an interaction term

We determined that I_crit is a more statistically valid performance indicator to evaluate CSS performance than the critical rainfall depth N_crit proposed in prior literature. Although the two proposed performance indicators appear about equally capable of predicting CSO occurrences (within the limitations of this study), I_crit is independent of fluctuations in wet weather from year to year. In addition, our multiple regression model predicts CSO volumes better than any bivariate correlation of CSO volume to an individual wet-weather variable, or other regression models considered.

Data limitations were encountered in this study, including the temporal and spatial resolution of rainfall data and metadata for CSOs. Rainfall data at 1-hour intervals were chosen because these data were the finest temporal scale data available over the entire 16-year period of analysis. While measurements at even finer time intervals could produce a multiple regression with greater ability to predict CSO volume, finer-scale temporal data were not available. Continuation flow data within each CSO outfall, connected impervious cover area per CSO outfall, amount of infiltration water, and storage volume per connected impervious area were also not available for the Cumberland CSS.

Cumberland is a small municipal system (treatment capacity 15 million gallons (∼56.8 million L) per day) with a relatively small, non-expanding urban area. As such, this assessment protocol may be less directly transferable to very large or very small systems that may experience different ranges of CSO discharge volumes and/or qualitatively different engineering issues, or that may be expanding or changing the amounts of impervious cover. We recommend replication of the model-building process in such systems to validate or modify the model and potentially controlling for changes in impervious cover, rather than extrapolating the regression of CSO discharge volumes well beyond the range of CSO discharge volumes in Cumberland.

Future work should consider what magnitude of change in the model's performance indicators – I_crit values, regression coefficients, and regression statistics – can be observed in association with documented past interventions to one or more applicable CSS(s) of interest. This approach (I_crit plus multiple regression) could help evaluate the effect of those interventions on reducing CSO incidents and volumes in a CSS. Alternatively, these performance indicators could be applied using future precipitation projections to identify how future changes in precipitation in a CSS watershed may impact CSS performance, assuming urbanization and the CSS remains relatively unchanged. Because this methodology is not dependent on the physical characteristics of an individual CSS itself, the modeling procedure could be scaled up to assess broader CSS performance across a group of systems or municipalities of similar size, to determine if CSSs in general are keeping up with wet-weather changes under climate change, and identify locales in need of assistance with doing so.

The authors gratefully acknowledge Citlalli Rojas Huerta and Larissa Chraim for data cleaning and initial analysis; Dr Eric Loken for regression assistance and helpful feedback; Dinah Voyles Pulver for assistance with data access and informative input; and Robert Smith, PE, Director of Engineering at the City of Cumberland, Maryland, for informative input and helpful feedback. Kirchhoff was partially supported by the National Science Foundation under Grant CMMI-1944664. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

All relevant data are available from an online repository or repositories (Bizer & Kirchhoff, 2022).

The authors declare there is no conflict.