Purpose and Scope

The purposes of this study are twofold: understand the drivers of historical water use in the Providence Water network (the service area of Providence Water and its wholesale customers) and forecast future water use in the same network. To meet these goals, a model was developed for quarterly water use, expressed as daily per capita rates of domestic, commercial, and industrial water use within each of the 8–10 PWSs connected to the Providence Water network. The model was developed for the span of January 2005–December 2021. This model was then used to make forecasts for each quarter in 2030 and 2040 under multiple simulated future scenarios. Quarters in this study run from January–March (Q1), April–June (Q2), July–September (Q3), and October–December (Q4). The scenarios simulated in this study included 12 combinations of 2 climate scenarios (low and high warming), 3 population scenarios (low, medium, and high growth rates), and 2 economic scenarios (low and high growth rates).

PWS Service Area Boundaries

In the following methods, a Geographic Information System (GIS) analysis was used to aggregate data for each PWS. The service area boundaries of each PWS were obtained from the Rhode Island Geographic Information System (RIGIS; Rhode Island Geographic Information System, 2022). For the service area boundaries of the two PWSs that later merged with Providence Water, an older version of the service area boundaries shapefile was obtained from the Rhode Island Water Resources Board (RIWRB). The service area boundary of Providence Water prior to these two mergers was determined by spatially subtracting out the service areas of the two merged entities.

Calculating Daily Per Capita Water-Use Rates

Early modeling work during this study identified population as the most influential driver of water use, and therefore water use was modeled as daily per capita rates. Modeling water use as per capita rates instead of as total withdrawals confers several advantages, including improved robustness to statistical methods and easier outlier detection (National Research Council, 2002). Historical data of daily per capita water-use rates (in gallons per capita per day [gal/cap/d]) were calculated from water sale volumes (in gallons) to customers for the three categories of water use previously described, from population served estimates for each PWS in the Providence Water network (in persons), and from the length of each quarter (in days). Average daily per capita rates were calculated for each quarter, year, PWS, and category of water use for which data were available between 2005 and 2021.

Predictor Variable Data

The meteorological predictor variables included daily maximum temperatures and daily precipitation values; these were obtained from the Global Historical Climatology Network (GHCN-DAILY; Menne and others, 2012, 2023). Temperature and precipitation data for the years 2000–22 within 25 kilometers of the study area (fig. 1) were filtered to remove data points with data quality flags. For each PWS, daily timeseries of both metrics were generated as a distance-weighted average of all values for stations within 25 kilometers of the centroid of each PWS service boundary. Daily temperature and precipitation data were then aggregated to quarterly values of mean maximum temperatures (in degrees Celsius) and total precipitation (in millimeters) for the years 2005–21.

Demographic predictor variables included population density and median age. Population density of each PWS was calculated for each year between 2005 and 2021 using the population served estimates and the PWS service areas. Median age was obtained from IPUMS at block-group level for years 2000–21 (Manson and others, 2023). Data gaps for each block group were filled in the way previously described in the “Population Served Estimates” section. If no values were available for a census block group, gaps were filled with the mean value of the predictor variable. Block group median age values were aggregated to PWS service area by weighted averaging, using the population of each block group as the weights.

Housing-related predictor variables included in this study were housing unit age, housing unit size, the percentage of units with complete kitchen facilities, and the fraction of housing structures that were single-family homes. The number of housing units, the median build year of housing structures, the median number of rooms per housing unit, and the number of housing units with complete kitchen facilities (including a sink with a faucet, a stove or range, and a refrigerator) were obtained for each block group from IPUMS, and gaps were addressed as previously described (Manson and others, 2023). The percentage of housing units with complete kitchen facilities was calculated for all block groups, and then the median build year, median number of rooms, and percentage of housing units with complete kitchens were aggregated to PWS service areas using weighted averaging. The number of housing units in each block group was used as the weights. The ratio of single-family housing structures was calculated for each PWS using the “E-911 Sites” data layer from RIGIS (Rhode Island Geographic Information System, 2021). This dataset is an emergency response tool that contains locations, descriptions, and classifications of all structures in Rhode Island as of 2021. Structure locations were assigned to PWS service areas using spatial joins, and the ratio of single-family structures (E-911 codes R1 and R3) to all residential structures (E-911 codes R1–R4) were computed for all PWSs. These ratios were considered static throughout the timeframe of this study.

Economic and labor predictor variables in this study included median household income, gross domestic product (GDP), number of private establishments, the coincident economic activity index (CEAI), the unemployment rate, and the percentage of the labor force in nine different sectors of the economy. Median household income data came from IPUMS at block group resolution (Manson and others, 2023), gaps were filled as previously described, and values were aggregated to the PWSs by weighted averaging using the number of housing units per block group as weights. Income data were adjusted for inflation to 2020 U.S. dollars. Inflation data were obtained from the World Bank application programming interface (API) through the “priceR” R package (Condylios, 2022; World Bank, 2023). Monthly data of the size of the labor force and the number of people employed for each county in Rhode Island were obtained from Federal Reserve Economic Data (FRED; Federal Reserve Bank of St. Louis, 2023) using the API service, gap-filled in the same way as the data from IPUMS and averaged by quarter. Quarterly unemployment rates for each county were calculated from these values. Unemployment rates for PWS service areas were the weighted average of the county values, weighted by the spatial overlap of each PWS and each county. Monthly statewide data on the number of employees in total nonfarming industries, construction, education and health services, financial activities, government, information, leisure and hospitality, manufacturing, mining and logging, and professional and business services were also obtained from FRED API (Federal Reserve Bank of St. Louis, 2023). They were averaged by quarter and converted to quarterly percentages of the labor force using the ratio of employees in each sector to the total number of employees in nonfarming industries. Annual GDP data and quarterly data of the total number of private establishments for each county in Rhode Island were obtained from FRED API (Federal Reserve Bank of St. Louis, 2023), gaps were filled as previously described, and data were normalized by the number of people employed in each county. Gross domestic product values were adjusted for inflation to 2020 USD (Condylios, 2022; World Bank, 2023). The coincident economic activity index (CEAI) is a composite index that includes the number of nonfarm employees, the unemployment rate, the average hours worked in manufacturing, and wages and salaries (Federal Reserve Bank of St. Louis, 2023). Monthly statewide CEAI values were obtained from FRED API and averaged for each quarter.

Model Development

A cubist regression machine-learning model was developed using the derived daily per capita rates and the corresponding predictor variable data (Quinlan, 1992; Eng and Wolock, 2022; RuleQuest Research, 2022). The cubist regression model is a type of decision-tree-based model that creates a list of rules about when to apply various linear regression equations to make model estimates of the data. The use of linear regression models allows for more nuanced model estimates under scenarios for which predictor variable data are outside the bounds of training data. Because this study used a model trained on historical data to forecast future scenarios in which conditions differ from those seen historically, this was an important distinction of cubist regression. Other common decision-tree models were considered, including Boosted Regression Trees (Friedman, 2001, 2002) and Extreme Gradient Boosting (Chen and Guestrin, 2016), but these two popular model structures split data into average values instead of linear regressions. As a result, forecasts made for the different future scenarios by these models were similar despite the differences between the future scenarios.

To prevent the model from overfitting the data (meaning that the model describes random fluctuations in the data instead of real signals and therefore cannot be generalized), the model was trained on 80 percent of the data (the training data), then had its predictive ability tested using the remaining 20 percent (the testing data). The daily per capita rates and the associated predictor variable values for each data point were split randomly into training and testing subsets. The modeling steps of selecting predictor variables to include in the model and tuning the model hyperparameters were performed using the training data. Once the model is tuned, the final version of the model is trained with the training data and tested on the testing data.

Cubist regression models have two hyperparameters called “committees” and “neighbors” (Kuhn and Johnson, 2013). Hyperparameters are variable parameters that determine how the model is structured. A committee is a way of boosting the model by building multiple models in sequence. In other words, the input data for each model are adjusted based on the residuals of the previous model, and the final model estimates are made by averaging over the committees. Specifically, for any mth committee, the model uses a pseudo-response variable (in this case, an adjusted daily per capita rate) of the following form:

where The number of committees determines the level of boosting applied to the model—more committees generally tend to produce better model estimates, but they can also quickly lead to overfitting the data.

When cubist regression models generate estimates, the estimates can be adjusted using the values of a certain number of neighbors, or similar data points from the training data (Kuhn and Johnson, 2013). This approach adjusts model estimates using the following formula:

where The inclusion of more neighbors in the estimation process leads to more smoothing of the model estimates. More details can be found in Kuhn and Johnson (2013).

The exact number of each hyperparameter used is determined during a model tuning process that tests combinations of hyperparameters. The combination found to produce a model that optimizes a given loss function is selected as the final “best” model. In this study, root mean square error (RMSE) was used as the loss function. Model tuning was done using tenfold cross validations (Kuhn and Johnson, 2013). The final combination of hyperparameters from the tenfold cross validation was determined using the “oneSE” approach, which avoids overfitting by selecting the simplest model within one standard error of the model with the most minimized loss function (Breiman and others, 1984).

Predictor variable selection was included in the modeling process to produce the most concise model that excludes predictor variables with negligible effect on daily per capita rates. The Recursive Feature Elimination (RFE) process from the “caret” package (Kuhn, 2008, 2019) was used for predictor variable selection. The process sequentially removes variables in reverse order of their importance, retrains the model, and determines performance using the remaining variables (Kuhn, 2019). Initially, all predictor variables were introduced, and each sequential elimination decreased the number of included predictor variables by five. Model tuning of hyperparameters was nested within the RFE process, with each round of predictor variable selection retuning a new model for that subset of predictor variables. The RFE process was performed within a tenfold cross-validation framework, and the subset of predictor variables retained in the final model was the smallest subset within one standard error of the model with the lowest RMSE (Breiman and others, 1984).

After tuning and recursive feature elimination, a final model was trained using the selected hyperparameters and predictor variables. Specifically, the model was trained with two committees, and predictions were made with three neighbors. Twenty-four predictors were introduced into the recursive feature elimination, and fifteen predictor variables were retained for the final model. After modeling, model estimates were post-processed to correct for systematic bias using the empirical distribution matching method (Belitz and Stackelberg, 2021). This method adjusts the values of the model estimates so that the distribution of the model estimates matches the distribution of the observed data.

Tuning and training were done using the R packages “caret” (Kuhn, 2008), “Cubist” (Kuhn and Quinlan, 2023), “foreach” (Microsoft Corporation and Weston, 2022b) and “doParallel” (Microsoft Corporation and Weston, 2022a), with data organization facilitated using the packages “dplyr” (Wickham and others, 2023a), “tidyr” (Wickham and others, 2023b), and “purrr” (Wickham and Henry, 2023). Investigation of the tuned model was performed using the “iml” package (Molnar and others, 2018). The historical data used to train and test the model, a script that trains the cubist regression model and generates model estimates, and the historical model estimates are available in Chamberlin (2024).

Forecasts

The cubist regression model was used to forecast future water use from simulated predictor variable data. Forecasts were made for each quarter of years 2030 and 2040 (a total of eight time periods) under 12 unique combinations of two climate scenarios (low and high warming), three demographic scenarios (low, medium, and high rates of population growth), and two economic scenarios (low and high rates of economic growth). To approximate uncertainty in forecasts, a Monte Carlo approach was taken in which 100 iterations of each scenario for each quarter were simulated. What follows is a description of how each predictor variable used for the forecasts was simulated under the various scenarios.

The fraction of single-family housing out of total housing, the percentage of housing units with complete kitchen facilities, and the median number of rooms per housing unit were kept the same as they were in the historical data for all future scenarios. Values for these prediction variables were simulated by randomly sampling the historical values of the prediction variable in the appropriate PWS. Median age of the population was simulated to increase at a similar rate of approximately 0.2 years per year under all the climate, demographic, and economic scenarios. This level of growth was selected based on the average rate of increase in the Nation’s median age between 2021 and 2022, after viewing historical data from the U.S. Census Bureau to confirm that this rate was consistent with historical growth (Manson and others, 2023; U.S. Census Bureau, 202346). Values for this prediction variable were simulated using the value for each PWS in 2020, and by adding to that a value sampled from a normal distribution with a mean of 2 and standard deviation of 1 for each decade. Labor forecasts similarly did not vary by scenario, but did vary between 2030 and 2040, and were based on projections made by the Rhode Island Department of Labor (Rhode Island Department of Labor and Training, 2023).

Mean daily maximum temperature for each quarter was varied by climate scenario to assess the effect of climate change on water use but was not varied between 2030 and 2040. This decision was due to the great amount of overlap in climate models between 2030 and 2040 for the northeastern United States (Hayhoe and others, 2007; Sun and others, 2015; Runkle and others, 2022). Mean daily maximum temperature values were simulated by sampling a normal distribution described by the mean and standard deviation of the historical mean daily maximum temperature values for each PWS and quarter, then adjusted by increasing the mean and standard deviations of the distributions by the values given in table 2. Adjustments were chosen using the work of Hayhoe and others (2007). In general, simulated temperatures are warmer under the high warming scenario: however, the variations around the means are large, and there is considerable overlap between the two scenarios.

Three population growth scenarios were simulated: high, medium, and low rates of growth. Populations served for each PWS were estimated for 2030 and 2040 using the 2020 population and simulated growth rates. Simulated growth rates were based on an analysis of historical growth rates in municipalities within the Providence Water network. Historical populations for 21 municipalities between 1970 and 2020 were obtained from IPUMS (Manson and others, 2023). For each of the 21 municipalities, five 10-year population change rates and four 20-year population change rates were calculated, producing a total of 105 historical 10-year population change rates, and 84 historical 20-year population change rates. The highest one-third of population change rates was randomly sampled to simulate population change in the high growth scenario, the middle one-third for the medium growth, and the lowest one-third for the low growth scenarios.¹ The distributions used are summarized in tables 3 and 4.

For each scenario, the population density of each PWS was calculated from the simulated population and the service area boundary. Service area boundaries were forecasted to remain constant. High rates of population growth were simulated to be accompanied by increased housing construction to accommodate larger populations, though the historical fraction of single housing out of total housing was kept constant. Median build year was increased from the 2020 value by a value sampled from a normal distribution with mean of 1 and standard deviation of 1 for each decade for the low growth scenario, a distribution with mean 3 and standard deviation of 1 for the medium growth scenario, and a distribution with mean 5 and standard deviation of 1 for the high growth scenario (tables 3, 4). These values were selected on the basis of a preliminary analysis of historical census data from IPUMS (Manson and others, 2023). Population growth scenarios were distinct from each other, with higher rates of population growth associated with higher total PWS populations served, greater average population densities, and newer average median build dates for housing structures (tables 3, 4).

Two economic scenarios were considered: high and low rates of growth. To calculate the future values of the unemployment rate, median household income, GDP per employee, number of establishments per employee, and CEAI, base values for all PWSs were adjusted using simulated change rates derived from the 2005–21 data used for model development. Ten-year rates of change for each predictor variable were calculated for each PWS. This produced 63 estimates of each predictor variable. For all predictor variables except unemployment, the higher half of estimates were sampled to simulate the high growth scenario and the lower half of estimates were sampled to simulate the low growth scenario. To simulate rates of unemployment, the higher half of estimates was sampled to simulate the low growth scenario, and the lower half of estimates was sampled to simulate the high growth scenario. These distributions are summarized in table 5. To simulate 2030 predictor variable data, distributions were sampled and applied to the 2020 base values once. To simulate 2040 predictor variable data, distributions were sampled and applied to the 2020 base values twice. Of note, is that the low growth scenario allows for simulating a prolonged recession, where unemployment increases, and median household income, GDP, and the number of establishments all decrease. The simulated values of the five predictor variables varied in the economic scenarios are summarized in table 5.

Model Structure and Performance

After model tuning and RFE, a cubist regression model with two committees was trained on the data. Model estimates were made using three neighbors. Of the 24 predictor variables introduced to the model, 15 were retained after the RFE process:

Data from the nine remaining predictor variables did not significantly reduce model RMSE, nor increase model r², and so they were excluded. The variables excluded were the fraction of the labor force in construction, education and health services, financial activities, government, information, leisure and hospitality, mining and logging, and professional and business services, and total precipitation.

After bias correction, model estimates of daily per capita rates were correlated with an overall Pearson’s r² of 0.99 for training data and 0.94 for testing data. Root mean square error for training data was 2.8 gal/cap/d and for testing data was 6.6 gal/cap/d. Modelled values of daily per capita rates for domestic, commercial, and industrial water-use fit the training data with Pearson’s r² values of 0.96, 0.97, and 0.80, respectively, and Pearson’s r² values of 0.77, 0.82, and 0.63, respectively, for the testing data. Root mean square error was 4.3, 1.9, and 1.1 gal/cap/d for training data for domestic, commercial, and industrial water use, respectively, and was 10.1, 4.6, and 1.3 gal/cap/d, respectively, for testing data. Model predictions are shown in figure 2.

Modeled total Providence Water network quarterly water use estimates from the cubist regression model included estimates of both reported and unreported water. Modeled total quarterly use between 2005 and 2021 had a significant negative trend of −0.96 Mgal/d per year (95-percent confidence interval of −0.79 to −1.15; seasonal Mann Kendall Trend test, p<0.001; fig. 3; Helsel and others, 2020). Estimated total Providence Water network use ranged from a minimum of 36 Mgal/d in Q1 2021 to a maximum of 73 Mgal/d in Q3 2005 (fig. 3). Estimated total Providence Water network water use should be higher than reported water use owing to the inclusion of unreported water use. It should also be less than total withdrawals from the Scituate Reservoir reported by Providence Water owing to water losses through leaks, water uses from the “Other (Total)” categories of the annual reports, water used to fight fires, and so on. Estimated total Providence Water network water use falls between these two values for 63 of the 68 quarters (or 93 percent of the period) examined in this study (fig. 3). Estimated water use fell below reported water use in Q4 2012, Q4 2015, and Q4 2020, and it rose above withdrawals from the Scituate Reservoir in Q4 2006 and Q1 2014. Modeled total Providence Water network quarterly use was significantly correlated with both reported use (p<0.001, Pearson’s r²=0.70) and Providence Water network withdrawals from the Scituate Reservoir (p<0.001, Pearson’s r²=0.46).

Shapley Additive exPlanations (SHAP) values show the relative effect of predictor variables within the cubist regression model for each model estimate (Lundberg and Lee, 2017). SHAP values that are greater than zero increase the estimate, and values less than zero decrease the estimate. The predictor variable with the highest average absolute SHAP values was “category of water use” (fig. 4). The designation of “domestic” was associated with higher water-use estimates, whereas designations of “industrial” or “commercial” were associated with lower water-use estimates (fig. 5). Quarters 1 and 2 are generally associated with lower water use than quarters 3 and 4, and the difference is much greater for domestic water use than for commercial or industrial (fig. 5). Temperatures above 20 degrees Celsius were associated with greater domestic water use (fig. 6). Higher population densities were associated with higher commercial water use and lower domestic water use (fig. 7). Housing units with more rooms, and a greater median year structure built (in other words, more recently built housing) were associated with higher commercial water use and lower domestic water use (fig. 8). The portion of total housing that was single-family homes did not have a monotonic effect (fig. 8). The relative variation in SHAP values between water-use categories was generally less for economic predictor variables than it was for demographic or housing-related variables. An exception is that higher rates of employment in the manufacturing sector were associated with greater increases in domestic water use than in commercial or industrial water use (fig. 9).

Future Forecasts

Forecasts made by the cubist regression model indicated that daily per capita rates are expected to decrease more in the scenarios with high population growth (fig. 10). However, total Providence Water network water use in Mgal/d is forecasted to be less sensitive to population growth than daily per capita water-use rates. Both daily per capita rates of water use and total Providence Water network water use are forecasted to decrease under higher rates of economic growth (fig. 11). Effects of climate change were minimal over the 2030 and 2040 prediction horizon based on the scenarios tested (fig. 12). Though the median water use forecasted is estimated to be slightly higher under higher warming scenarios, the range of forecasted water use under each climate scenario overlap substantially.

Across all scenarios, annual average water use was forecast to decrease in 2030 and further decrease in 2040 (fig. 13). The annual forecasted daily per capita rates and Providence Water networkwide water use are forecasted to be significantly below the levels of water use seen in 2005–21 (Kruskal-Wallis test, p<0.001 for both groups). Of all the forecasts, only 17 percent of forecasted daily per capita rates and 21 percent of forecasted annual Providence Water networkwide uses were equal to, or greater than, those seen in 2005–21. Of all forecasts, 83 percent of daily per capita rates and 79 percent of Providence Water networkwide water-use forecasts were below values from 2005 to 2021.

Annual median forecasts of daily per capita rates varied from a low of 42.4 gal/cap/d in 2040 under high population and economic growth scenarios and high climate warming to a high of 67.6 gal/cap/d in 2030 under low economic and population growth scenarios with high climate warming (table 6). The lowest median annual forecasted Providence Water networkwide use was 35.0 Mgal/d in 2040 under medium population growth, high economic growth, and low climate warming. The highest median annual forecasted Providence Water networkwide water use was 43.5 Mgal/d in 2030 under low population and economic growth, and high climate warming (table 6). Median forecasts for all scenarios were less than model estimates for 2005–21. The 95-percent confidence intervals of forecasted scenarios had some overlap with the historical 95-percent confidence intervals for all but the daily per capita water-use rate forecasts in the high population growth, high economic growth, and high climate change scenario in 2030. For the 2040 forecasts, however, only low population growth scenarios had overlapping 95-percent confidence intervals with the historical estimates for the daily per capita water-use rates (table 6). In 2040, all the high population growth rate scenarios had overlapping 95-percent confidence intervals for the total Providence Water network water use.

Similar patterns of forecasts were observed at the quarterly timescale, where forecasts fell below the historical ranges of quarterly water use (fig. 14). The differences between forecasted water-use rates and simulated water-use rates for the years 2005–21 were comparable in magnitude to the differences between quarters for each forecasted year. For example, daily per capita rates in Q3 2040 were predicted to be similar to the rates in Q1 each year from 2005 to 2021 (fig. 14). Forecasts for Q2 and Q4 were similar to each other in the forecasts for 2030 and 2040, whereas the Q3 forecasts were higher, and Q1 forecasts were lower. Quarterly forecasts are summarized in table 7.

Drivers of Historical Water Use

The seasonality of water use, especially of domestic water use, is seen in many water-use studies, with more outdoor water use during the summer months (Mills and others, 2014; Ahmed and others, 2020). The SHAP values for the predictor variables average daily maximum temperature and calendar year quarter agree with this pattern: water use at temperatures greater than 20 degrees Celsius was greater than water use at lower temperatures, and the greatest increase was for domestic water use (fig. 6). The importance of calendar year quarter as a predictor variable in the model, with the highest water use in summer and fall, also supports greater seasonality for domestic water use than for commercial or industrial water use in this study (fig. 5).

Outdoor water use would reasonably be expected to be higher for single-family homes than multi-family homes, leading to greater overall daily per capita water-use rates in single-family homes (Stoker and Rothfeder, 2014; Villarin and Rodriguez-Galiano, 2019). In this study, areas with comparatively high proportions of single-family homes did have higher residential water use; however, so did areas with comparatively low proportions (fig. 8). The proportion of single-family homes was the only variable that was fixed for each PWS and did not vary across quarters or years; therefore, it is possible that the model used this predictor variable to account for differences between PWS not explicitly accounted for in this study (for example, differences in reporting procedures or differences in how service connections are categorized). The increase in water use at lower proportions of single-family housing may be due to some other factor not accounted for in this study. This study found a decrease in daily per capita domestic water use with increasing population density (fig. 7), which was consistent with the pattern of decreased per capita water use for larger cities seen in the United States (Mahjabin and others, 2018).

A previous study of water use in Providence, Rhode Island, done by Stagnitta and Medalie (2023), used linear regression to model water use at a monthly timescale between 2014 and 2021. The results of this study were generally consistent with Stagnitta and Medalie (2023). Their study also modeled daily per capita rates, and modeled single-family domestic, multi-family domestic, commercial (defined slightly differently than this study), and industrial water use separately with unique linear regression models for each category of water use. For both studies, temperature (as average daily maximum temperature in this study and as the square of minimum temperature in Stagnitta and Medalie [2023]) had a positive relation with domestic and industrial water use, though temperature was not retained in the commercial water-use model of Stagnitta and Medalie (2023). Similarly, Stagnitta and Medalie (2023) found binary indicators of low and high water-use seasons (defined as December–May for low water use and July–September for high water use) to be influential in their study. These findings may correspond to the influence of the calendar year quarter predictor variable in this model. Both Stagnitta and Medalie (2023) and this current study found CEAI to be inversely related to water use. The effects of the coronavirus disease 2019 (COVID-19) pandemic were explicitly included in Stagnitta and Medalie (2023) and were found to be significant for the commercial and industrial models. In the current study, the effects of the COVID-19 pandemic were not explicitly included, but the model estimated water use in 2020 (with an r² of 0.98 on training data and an r² of 0.88 on testing data) using the provided input predictor variable data. These model results may indicate that the provided input data sufficiently represented the effects of the pandemic without having to include the pandemic explicitly.

Forecasts of Future Water Use

The model estimates for many future scenarios showed water-use values that were less than reported or estimated water-use values over the past two decades. Median water-use forecasts were generally higher for 2030 than 2040, and 2030 projections were mostly lower than the water use estimated for 2021. The daily per capita rates are forecast to decrease with increasing population (due to increasing population density and newer housing) however, increasing populations increase total water use (in that more people consume more water), and, as such, the relative balance of these two effects led to a nearly negligible effect on total Providence Water network water use according to the population scenarios (fig. 10). Historically, the increase in population over the past two decades has been associated with an overall decrease in water use, which is likely caused by increased water efficiency of domestic use (DeOreo and Mayer, 2012; DeOreo and others, 2016). The model did not directly account for increases in efficiency of water-using appliances; however, other predictor variables, such as the median build age of housing units, may be capable of representing this trend. Other factors that can improve water efficiency, such as controlled leakage in pipes, more effective water auditing, and consumer education campaigns (U.S. Environmental Protection Agency, 2016), were likewise not included in this model, because of a lack of representative data.

The historical decrease of total Providence Water network water use between 2005 and 2021 indicates that the decreases in daily per capita rates have more than compensated for increases in population. The forecasts from this model indicate that this balance might be shifting in the future, in that decreases in daily per capita rates may no longer outweigh increases in population, and total network water use may not continue decreasing at the same rate or might stabilize at a lower level of water use. Though nearly all simulated scenarios were forecasted to have lower water use than historical water use, 17 percent of forecasted daily per capita rates were equal to historical rates. Therefore, the forecasts leave open the possibility for an increase in future water use, though the likelihood of a future increase in water use is small.

A limitation of this study is that the cubist regression model assumes that the current associations between the data will continue into the future. For example, if the association between changes in housing stock and changes in water use are different in the future (for example, if daily per capita rates no longer decrease with increases in median age of houses), the model will not account for these shifts. This study attempted to account for future unknowns such as this by modeling multiple future scenarios. We recommend using the forecasts based on demographic, economic, and climate scenarios together to bound the expectations for future water use projections rather than using them in isolation as precise estimates for individual scenarios.