Monthly-to-Daily Water-Use Data Comparison

Prior to using these data to build models and obtain predictions, NJAW’s daily public supply data were first compared to monthly public supply water-use data from the state-owned NJWaTr database (New Jersey Department of Environmental Protection Division of Water Supply and Geoscience, 2022). The purpose of this comparison was to help validate both datasets and identify any inconsistencies that may affect the analysis and models.

Daily system delivery data acquired from NJAW (in units of million gallons per day) were temporally aggregated to monthly values and then compared to monthly DWSA system data from the NJWaTr database for the 15 NJAW DWSA systems for the years 2016–20 (figs. 3, 4, and 5). Of the 15 DWSA systems, nine showed a close match to the monthly NJWaTr public supply data (table 1). In some cases, the water-use data from the two datasets matched almost perfectly, so that they were often indistinguishable when plotted, as is the case with Ocean City (fig. 3B), Strathmere (fig. 3C), Washington-Oxford (2016–19; fig. 4E), Harrison (fig. 5C), Bridgeport (fig. 5D), and Logan (fig. 5E) DWSA systems.

Three of the 15 DWSA systems (referred to in this study as Barrier Islands [fig. 3D], Burlington [fig. 5A], and Camden [fig. 5B]) are portions of larger DWSA systems for which complete system data were not provided by NJAW: the Coastal North (part of the Coastal Region [NJ134500]) and Western Division (part of the South Region [NJ0327001]) systems. For this study, they were treated as individual systems. Because these three systems represent only partial system data, the aggregated daily data values from NJAW did not match well with, and were much lower than, the corresponding NJWaTr monthly values for the Coastal North and Western Division systems with a mean percent difference of at least −134% or smaller. For the remaining three systems, Passaic, Little Falls, and Frenchtown, the data from NJAW and NJWaTr do not match well even though they were complete system data and not partial system datasets. In these three systems, it is likely that the aggregated daily data and the monthly NJWaTr are different because of reporting inconsistencies, and (or) the complex nature of systems that exist in the northern part of the state (Vince Monaco, NJAW, oral commun., June 2020).

Daily Public Supply Water-Use Data

Daily public supply water use in New Jersey varies by time of year and location. Generally, publicly supplied water use is greater in the summer months than in the winter months (NJDEP, 2017), although this was not always the case for this study. Across all 15 systems, the average water use per day in the summer months (June–August) was 116 percent greater than the winter months (December–February). From visual examination of the plotted daily public supply water-use data, it appeared that a few data points were noticeably outside the typical range of values for the Strathmere, Frenchtown, Burlington, Camden, and Bridgeport systems (figs. 6C, 7B, 8A, 8B, and 8D). For these 10–20 data points, plant engineers from NJAW were consulted. Where possible, they provided explanations for what occurred at those plants on some of those days of interest. These explanations included hydrant flushing, water main breaks, using water to fight fires, or water-quality testing of newly replaced mains (NJAW, written commun., 2021). Some further analysis was conducted to estimate the effect of these data points by comparing the model results based on the datasets with and without the data points of interest. In addition to these short-term anomalous data points, the Washington-Oxford system had a noticeable shift in water use starting around the beginning of 2019 (fig. 7E). The average water use for this system from 2016 through 2019 was 1.10 million gallons per day (Mgal/d); thereafter, the average water use was 2.19 Mgal/d through 2020. This apparent shift in water use was possibly caused by changes in the DWSA system boundaries or changes in the buying and selling of water to other DWSA systems; however, the exact cause of the shift in water usage is unknown. Another noticeable long-term change in daily water use is observed in the Little Falls system about mid-way through 2018 (fig. 7C). The magnitude of volumes does not noticeably shift, rather the day-to-day variability appears to increase. This change may be due to differences in the way the data were reported. Finally, the daily dataset provided for the DWSA system of Frenchtown only included years from 2017 through 2020 and, as a result, the analyses for this system are based on 4 years of data. These noticeable shifts in data for the Washington-Oxford and Little Falls systems and missing data for the Frenchtown system occur over a longer timeframe in contrast to systems with anomalous data which occur over only 1 or 2 days. As a result, these longer-term shifts or data gaps would likely have a larger effect on any analysis and modeling performed as compared to those systems with shorter-term anomalous data because there are more anomalous data points to influence the model. However, for this study, the same analysis was carried out among all systems regardless of the apparent long-term shifts in datasets and missing data.

The 15 systems were grouped into two categories based on general, observable patterns in the daily data. Nine systems showed a strong seasonal signal throughout the 5 years of data, where water use was higher in summer months compared to winter months (figs. 6A, 6B, 6C, 6D, 7A, 8A, 8C, 8D, and 8E). The remaining six systems showed little to no seasonal pattern; instead, daily water use was mostly constant throughout each year (figs. 7B, 7C, 7D, 7E, 8B, and 8F). Among the nine DWSA systems with a seasonal signal, or seasonal systems, the average water use was 185-percent higher in summer months (June–August) compared to winter months (December–February), whereas the six DWSA systems with little to no seasonal signal, or non-seasonal systems, had an average increase in water use of 13 percent in summer months over winter months.

The initial groupings were confirmed by identifying the periodicity of each system’s data through periodograms. Periodograms quantify the relative importance of different frequencies in a dataset based on a scaled-Fourier transform of the time series data (Bloomfield, 2000). Peaks or local maxima in a periodogram indicate dominant frequencies in the data. Identifying the frequency at which the maximum spectral density occurs then allows for the corresponding period to be determined by computing the inverse of the frequency of interest, where the spectral density is the relative strength of the frequencies within the given time series signal (Kendall, 1946; Iyer and Chowdhury, 2009). Raw periodograms are rough estimates of the true spectral density (frequency domain representation of time series data) and are therefore subject to fluctuations and noise. Applying smoothing can improve the stability of (or the extent to which unimportant details are removed from) the raw periodogram (Bloomfield, 2000). The modified Daniell filter is one of the most common methods for smoothing, and in essence, constitutes a weighted moving average window where the endpoints receive less weight than the interior points. The modified Daniell filter can also be applied successively for more extensive smoothing (Bloomfield, 2000).

By generating a smoothed periodogram for each of the 15 systems, the seasonal and non-seasonal groupings were identified quantitatively (fig. 9). Systems with periodograms that had a maximum spectral density of approximately 365.2 days were categorized as seasonal systems (figs. 9A, 9C, 9D, 9E, 9F, 9G, 9I, and 9K), whereas systems with periodograms that had a maximum spectral density corresponding to a different number of days (less than or greater than 365.2) were categorized as non-seasonal (figs. 9B, 9H, 9J, 9L, 9M, 9N, and 9O). The groupings from the periodogram analysis confirmed the preliminary visual-based groupings, identifying the same nine systems as seasonal and the same six systems as non-seasonal. Burlington, Cape May Courthouse, Ocean City, Strathmere, Passaic, Harrison, Bridgeport, Logan, and Barrier Islands were confirmed as the nine seasonal systems. Camden, Frenchtown, Little Falls, Penns Grove, Belvidere, and Washington-Oxford were confirmed as the six non-seasonal systems. These seasonal and non-seasonal groupings were the basis for developing two regression equations to account for this key difference in the systems: one model equation form was developed for seasonal systems, and one model equation form was developed for non-seasonal systems.

Datasets and Methods

The datasets used in characterizing the DWSA systems in New Jersey included population served by DWSA for 2021 from the EPA’s Safe Drinking Water Information System (SDWIS; U.S. Environmental Protection Agency, 2021), the 2010 population estimates by census block group from the U.S. Census Bureau (U.S. Census Bureau, 2010), the 2017 New Jersey geographic boundaries for all DWSA systems (NJDEP Bureau of GIS, 2017b), the state of New Jersey’s 2019 tax parcel and property value data (New Jersey Office of Information Technology Office of Geographic Information System, 2019), the 2015 land-use and land cover data (NJDEP Bureau of GIS, 2015), the 2019 median household income estimates by census block group (U.S. Census Bureau, 2019), the landscape regions of New Jersey (NJDEP Bureau of GIS, 2017a), and NJWaTr monthly water-use data (New Jersey Department of Environmental Protection Division of Water Supply and Geoscience, 2022).

All datasets were clipped and summarized at the DWSA system level. The DWSA boundary coverage shapefile, obtained in 2020, represented the 2017 DWSA system boundaries and contained 589 unique systems (NJDEP Bureau of GIS, 2017b). The population-served data from SDWIS were available at the DWSA system level. Three of the DWSA systems from the NJAW daily datasets were provided as partial systems; to keep the analysis as consistent as possible, these partial systems were treated as individual systems. To help obtain population-served estimates for these partial systems, the 2010 Census Bureau population estimates were also used. The proportion of the total population in the partial DWSA system to the total population in the complete DWSA system was applied to scale SDWIS population-served data, resulting in estimated partial populations served for the three partial DWSA systems.

Median property values and median household incomes by DWSA system were calculated by aggregating land parcels and census block groups, respectively, to the DWSA system level. For these datasets, land parcels and census block groups were first intersected with the DWSA system boundaries. The property values dataset was filtered by property class, so only Class 2, or “residential property,” was considered for this analysis (New Jersey Register, 2018). Additionally, land parcels with a property value of zero were excluded. The median property values and household income were then calculated for each DWSA system. Because the data were filtered to include only residential properties, some of the 589 unique DWSA systems did not have appropriate tax data and thus did not have a median value. For these few systems, they were removed from the analysis as there were no data available. Median, rather than mean, values were used for the characterizations because the data were not normally distributed.

After obtaining a single summarized value for each socio-economic dataset (population served, household income, property value) per DWSA system, quartiles were calculated for each dataset based on all DWSA systems combined. Using quartiles allowed for the systems to be categorized into four groups based on the value of a given system. Subsequently, a given system was grouped into quartile 1 (Q1), quartile 2 (Q2), quartile 3 (Q3), or quartile 4 (Q4) according to the median property value for that particular DWSA system. This method of calculating quartiles to group the systems was used for the population-served values, median property values, and median household income values datasets (table 2).

The residential land-use category was the only category from the 2015 land-use data used for this analysis, which accounts for about 12.5 percent of total land-use area in New Jersey (NJDEP Bureau of GIS, 2015). The residential land-use category is further separated into residential density subcategories. The rural residential density land-use subcategory comprises 13.8 percent, the low residential density land-use subcategory comprises 18.5 percent, the medium residential density land-use subcategory comprises 48.6 percent, and the high residential density land-use subcategory comprises 19.1 percent (NJDEP Bureau of GIS, 2015). These data were first intersected with the DWSA system boundaries, then residential density percentages were calculated at the DWSA system level. Initially, the percentage of land area for each of the four residential density categories (rural, low, medium, and high) was calculated for each system, where the total area of land for each density subcategory was divided by the total residential land-use area (in acres). To designate a single value for each DWSA system, the residential density with the largest percentage was used to categorize each system. Because there was little observable difference in water-use values with respect to seasonality between DWSA systems with rural, low, or medium residential densities as the largest percentage, these three residential densities were combined and classified as “non-urban.” Systems with high residential density as the largest percentage were classified as “urban.” The distinction between urban and non-urban systems was used in this study because there was a difference in seasonality between systems with rural, low, or medium residential densities and systems with high residential densities.

Another factor that was used to classify the DWSA systems was whether a given system was coastal or not. In New Jersey, many of the coastal towns and areas receive a large population influx in the summer from tourists (Stirling, 2018). Obtaining estimates on changes in these populations over the seasons was difficult, so identifying whether a DWSA was geographically close to the Atlantic coast was used as a proxy for the change in population caused by tourism. To determine a DWSA system’s status as coastal or not, the NJDEP landscape regions of New Jersey (NJDEP Bureau of GIS, 2017a) was used in conjunction with the DWSA system boundaries (NJDEP Bureau of GIS, 2017b). The “Atlantic coastal” region definition from the dataset provides a reasonable coverage of coastal, tourist-destination towns in the state (fig. 10). The DWSA systems that intersected with the Atlantic coastal landscape region were identified and categorized as coastal systems.

Alternative methods to identify and categorize coastal systems were explored but not used because the categorizations of coastal systems included systems outside of the typical New Jersey summer tourist destinations, such as along the Delaware, Newark, and Raritan Bays. For example, the Coastal Area Facilities Review Act (CAFRA) boundaries for New Jersey dataset (NJDEP Bureau of GIS, 2005) created for planning and permitting purposes, represents coastal planning areas, but was not selected because it included DWSA systems along the Delaware Bay and along the Delaware River as well as some systems over 12 miles inland from any coastline. Another method considered was to apply a 5-mile buffer to the U.S. country boundary along the New Jersey coast. These alternative datasets were intersected with the DWSA boundaries to identify potential coastal systems. In both cases, there were DWSA systems included in regions not typically considered coastal towns or summer tourist destinations such as along the Raritan Bay (including Perth Amboy to Newark, New Jersey) and on the southern end, along the Delaware Bay to Salem County. The landscape regions of New Jersey dataset was used to assign coastal designations to the different systems so that systems that were not representative of the tourism effect in New Jersey and were not in the near vicinity of the Atlantic coast were not included. The landscape regions dataset best identified systems encompassing summer destinations in the state, that are affected by tourism and result in seasonal population changes, which impacts water-use patterns throughout the year.

Identifying Seasonality in Monthly Data

Once all DWSA systems were characterized into socio-economic (population served by DWSA, household income, property value) and geographical (coastal and non-coastal, urban and non-urban) groupings, monthly water-use data from NJWaTr were analyzed for seasonality (appendix 1). To identify whether a given system had a seasonal pattern, at least 24 consecutive months of water-use data from 2016 through 2020 were required. These years were the most recent 5 years of available data in NJWaTr at the time of analysis and most closely aligned with the daily data obtained from NJAW. Because of this requirement, 158 DWSA systems were excluded from this analysis, and 434 DWSA systems remained. Seasonality in monthly water use was determined using the same method described in the “Daily Public Supply Water-Use Data” section and is defined as higher usage during the summer months (May–September), as compared to lower usage during the cooler fall, winter, and spring months (October–April). The periodograms generated for Brigantine WD, Egg Harbor City WD, Norms Dale MHP, and Verona WD show how the strength of frequencies varied among the 434 systems (figs. 11 and 12). Systems with strong seasonality, like Brigantine WD and Egg Harbor City WD, show a peak frequency at 1/12, or a peak period of 12 months (figs. 11A, 11B, 12A, and 12B). It was determined that data could still visually appear seasonal with a dominant period of anywhere between 11 and 13 months, like that of Verona WD; therefore, any system with a peak in the periodogram corresponding to 11 through 13 months was classified as seasonal (figs. 11D and 12D).

In addition to using periodograms, the monthly water-use data for each of the DWSA systems were plotted and visually examined for a seasonal signal. There were 69 instances where the visual check was inconclusive so the periodogram results were used. There were 15 instances where the visual check and periodogram results disagreed. There were five cases where the visual inspection indicated that the systems were not seasonal, despite the periodogram results indicating they were seasonal. The same inspection found 10 systems were seasonal, despite the periodogram results indicating they were not seasonal. When the visual check was not immediately clear on whether a system was seasonal or not, the automated periodogram results were used. Only when a system’s water-use patterns were visually very clearly seasonal or not, was the visual check used in place of the periodogram as the periodogram can be sensitive to noise (one or two months with anomalous data) or large shifts in overall magnitudes in the data (a system’s usage noticeably decreases or increases over time). Another limitation of the periodogram is that it can only identify if there is any type of repeating, periodic pattern and does not distinguish between which type of periodic pattern. In this case, a specific seasonal pattern is the only repeating pattern of interest, where there is a substantial increase in water use in the summer months as compared to winter months. Occasionally, the results of the periodogram would indicate there is a yearly (or 12 month) repeating signal, but it would not follow the specific seasonal pattern described above. Because of these limitations, the visual inspection was used to verify or check the results.

Results and Patterns in DWSA Characterizations

One of the main purposes of the DWSA characterization was to identify any patterns in characteristics across all the systems, as they related to the degree of seasonality in water use. The motivation for focusing on seasonal water-use patterns came from the noticeable distinction in the daily data from the 15 NJAW systems highlighted in this study. Most of the 15 systems were placed into either a seasonal or non-seasonal group based on the daily water-use data provided. Therefore, the ability to determine if a system has a seasonal water-use pattern helps guide model development and water-use estimates for any DWSA system. All 15 NJAW DWSA systems were characterized based on the five datasets discussed above and were grouped into seasonal and non-seasonal categories based on the daily data provided from NJAW (table 3). The remaining 434 DWSA systems with water-use data in NJWaTr spanning the years 2016–20 were characterized in the same manner using monthly NJWaTr data to assess seasonality (appendix 1).

The DWSA characterization results were graphed by the five datasets broken up into quartiles (figs. 13, 14, 15) or qualitative categories (fig. 16). Generally, there are more DWSA systems that display seasonal water-use patterns for those systems with larger populations served, higher median household income, and higher median property values as compared to smaller populations served, lower median household income, and lower median property-value systems. This can be observed by comparing the percentage of seasonal and non-seasonal systems between the highest quartiles (Q4) and lowest quartiles (Q1; figs. 13, 14, 15). This pattern is most distinct in the population-served dataset. In Q4 for population served, 80 percent of systems showed seasonality, whereas in Q1, only 20 percent of systems showed seasonality, a decrease of 60 percent (fig. 13). This pattern is least observable in the household income dataset where the corresponding difference between Q4 (64 percent) and Q1 (50 percent) is a decrease of only 14 percent (fig. 14). Overall, the percentage of seasonal DWSA systems varied the least when grouped by household income, which suggests this variable may not be as influential as others when it comes to predicting seasonality in monthly water use for New Jersey.

When systems are grouped by coastal and non-coastal categories, it is evident that more coastal systems show seasonality than non-coastal systems (fig. 16B). This finding supports what is known about increased summer populations in coastal areas, where increased populations during summer months are associated with increases in publicly supplied water use for those same months. When systems are grouped by urban and non-urban residential density, urban systems often show non-seasonal water-use patterns.

Ninety-seven percent of all DWSA systems that were classified as coastal showed a seasonal signal. Seventy percent of all systems classified as non-urban showed a seasonal signal. Additionally, 80 percent of all systems in Q4 for population served showed a seasonal signal (fig. 17). These results indicate that a DWSA system shows seasonal patterns of water use if it is classified as coastal, non-urban, and (or) it serves a large population (Q4).

The systems were then grouped into more specific categories based on two characteristics rather than just one. Grouping the systems in two categories provided insight on whether certain DWSA characteristics had a stronger or weaker influence on seasonal water-use patterns. Combinations based on two of the three characteristics highlighted in figure 17 were used to create a total of nine distinct groups (figs. 18, 19, and 20). To align the population-served category with the qualitative categories (coastal or non-coastal, urban or non-urban), the population-served category was split into two groups: the top 50 percent (Q3 and Q4; denoted as Q3+Q4) and the bottom 50 percent (Q1 and Q2; denoted as Q1+Q2).

When the coastal systems are further grouped by residential density and by population served, they all show a high percentage of seasonality. When categorized in this manner, using these four categories (urban, non-urban, top 50 percent population [Q3+Q4], and bottom 50 percent [Q1+Q2] population), 82 percent or more systems displayed seasonal water-use behavior as seen in figure 18. This suggests that a system’s proximity to the coast is more influential in predicting its seasonality compared to the other factors considered here, namely residential land-use type and population size.

Systems that serve the top 50 percent populations and are categorized as coastal show the strength of the coastal factor as a determinant in seasonality (fig. 19). Seasonality was observed in 100 percent of coastal systems that serve top 50 percent populations, compared to 69 percent of non-coastal systems that serve top 50 percent populations (fig. 19). When separated by urban and non-urban residential densities, the percentage of seasonal systems that serve top 50 percent populations were 68 percent and 78 percent, respectively. The difference between these seasonal systems also indicates that the residential density is influential in predicting whether a system exhibits seasonality, where urban systems serving large (or top 50 percent) populations are less likely to exhibit seasonality than non-urban systems serving large (or top 50 percent) populations (fig. 19).

Separating urban systems by coastal classification shows that the coastal factor is likely a stronger determinant for a system’s seasonality than the residential density factor. Across all urban systems, 42 percent showed a seasonal signal (fig. 16A). But, when separated by coastal classification, 96 percent of urban, coastal systems showed a seasonal signal and 29 percent of urban, non-coastal systems showed a seasonal signal (fig. 20). The large difference in percentages when grouped further by coastal and non-coastal systems indicates a system’s proximity to the coast is likely a more influential determinant for its seasonality than its degree of urban residential density. If residential land use was more influential, the percentage of seasonal urban, coastal systems and seasonal urban, non-coastal systems may be more similar, or the percentage of seasonal urban, coastal systems may be much lower than the overall coastal system percentage (97-percent seasonal; fig. 17). The similarity in the percentages of seasonal urban, coastal systems and seasonal non-urban, coastal systems (98 percent; fig. 20) further indicates that a system’s coastal classification is more influential than a system’s residential density classification. Increased water use in summer months in systems along the coast, or in closer proximity to the coast, may be due to the increased number of people in those areas as a result of the influx of summer residents and tourists (Stirling, 2018). Additionally, higher water usage in summer months may be further explained by lawn irrigation and other outdoor activities such as car washing, flower and shrub watering, and use of sprinklers or pools (Dieter and others, 2018).

Similarly, there is an observable difference in seasonality percentages between urban systems that serve top 50 percent populations (68-percent seasonal) and urban systems that serve bottom 50 percent populations (19-percent seasonal; fig. 21). The percentage of systems that are seasonal between urban systems that serve top 50 percent populations (68 percent) and non-urban systems that serve top 50 percent populations (78 percent) is comparable and may indicate that the size of a population served may be more influential than the system’s residential density classification (fig. 21).

In summary, the degree to which a system’s water-use patterns show seasonality over a 12-month period, as defined as higher usage during the summer months (May–September), as compared to cooler fall, winter, and spring months (October–April), seems to be highly dependent on proximity to the coast and size of population served. If a system is classified as urban, it is more likely to display seasonality if the system serves a large population or is near the coast. Furthermore, if a system is classified as non-urban, it is likely to show seasonality, and even more so, if it serves a large population and (or) is near the coast.

Datasets Incorporated

Daily water use varies day to day and is influenced by many factors. Weather variables, such as temperature, precipitation, and evapotranspiration, are often considered influential factors affecting daily changes in water use; numerous studies have evaluated their relation to daily and monthly water-use estimation (Maidment and Miaou, 1986; Zhou and others, 2000; Eslamian and others, 2016; Opalinski and others, 2019; Ahmed and others, 2020). In addition to weather-related variables, other temporal variables have often been included in daily water-use estimation studies, such as season, day of the week, or whether the day falls on a weekday or the weekend, and holidays (Wong and others, 2010; Eslamian and others, 2016). There are also factors that affect water use on longer time scales, such as year to year. Socio-economic variables such as household income, price of water, and population change, are often considered when studying annual, or other longer-term changes in water use (National Research Council, 2002; Eslamian and others, 2016). Variables affecting long-term water use were not considered in the model development because the primary purpose of this study was to estimate daily water use. However, some of the variables considered influential in long-term water-use trends were used to describe and categorize all the DWSA systems in New Jersey (see previous section titled “Drinking Water Service Area System Characterizations”).

To estimate daily public supply water use in New Jersey, multi-variable linear regression models, using ordinary least squares regression methods, were developed for each of the 15 DWSA systems. Daily maximum temperature, daily precipitation total, number of days since significant precipitation (defined as greater than 0.15 inches), reference evapotranspiration¹ (ET), season, and day of the week were all considered in the model development. These initial variables were chosen based on review of previous studies and based on datasets that were readily accessible (Wong and others, 2010; Eslamian and others, 2016; Ahmed and others, 2020). Daily maximum temperature in degrees Celsius, daily precipitation in millimeters, and daily reference evapotranspiration in millimeters per day were obtained and downloaded from gridMET, a gridded dataset (with an approximate 4-kilometer spatial resolution) of meteorological data available through the Climatology Lab (Abatzoglou, undated), using climateR (version 0.1.0; Johnson, 2021). The dataset spans the contiguous United States from 1979 onwards, updated daily (Abatzoglou, 2011). The number of days since significant precipitation data was calculated from the daily precipitation totals to tally days since significant precipitation, using data downloaded from gridMet (Abatzoglou, 2011; undated). The season and day-of-the-week datasets were based on the time period of the daily public supply dataset (2016–20). The four seasons were defined as winter (January–March), spring (April–June), summer (July–September), and fall (October–December).

The DWSA boundary dataset for New Jersey (NJDEP Bureau of GIS, 2017b) was used to compute the mean daily maximum temperature and mean precipitation totals by DWSA in degrees Fahrenheit and inches, respectively. The volume threshold for the days since precipitation metric was set at 0.15 inches per day. Whether the day fell during the week or on the weekend was information used as a predictor variable, instead of the specific day of the week, because daily public supply data are known to be influenced by weekday or weekend (Eslamian and others, 2016). This factor is heretofore referred to as the weekday-weekend effect.

All analyses were performed using R statistical software (version 3.6.1; R Core Team, 2019). All potential predictor variables were first examined for multicollinearity, a statistical concept in which two or more predictor variables are highly correlated with one another (Helsel and others, 2020). Multicollinearity was estimated using the variance inflation factor (VIF) metric. Daily reference ET rates and daily maximum temperature showed high multicollinearity. For linear regression methods, it is assumed that all predictor variables are independent of each other and not highly correlated. As daily maximum temperature is more easily measured (and thus, likely more accurate) compared to ET, and because temperature is one of the most common predictor variables used in water-use estimation studies, ET was removed from consideration in building the regression models and only temperature was retained. This helped avoid the issue of ET and temperature being highly correlated to each other and not independent variables.

The data were also tested for normality using the Shapiro-Wilk test (Helsel and others, 2020) as well as visual examination of the distribution of the data. Because it was found that the data were not normally distributed, correlations between each predictor variable and the response variable (daily public supply water use) were estimated using the Spearman’s correlation coefficient where a p<0.05 indicates a 95-percent confidence level of correlation between the two variables (Helsel and others, 2020). Spearman’s correlation coefficient is nonparametric and therefore does not rely on any assumptions about the distribution of the data. The results of this correlation analysis indicated that season was not significantly correlated with daily public supply water use; thus, season was removed as a predictor variable and was not considered further in this study.

All nine NJAW DWSA systems categorized as displaying seasonal water-use patterns showed a significant correlation (p<0.001) between daily public supply volumes and daily maximum temperature (table 4). Within the same group, eight out of the nine seasonal systems showed significant correlation (p<0.05) between daily public supply volumes and the weekday-weekend effect, six systems showed significant correlation between daily public supply volumes and days since significant precipitation (p<0.01), and two systems showed significant correlation between daily public supply volumes and daily precipitation totals (p<0.05; table 4). Within the six NJAW DWSA systems categorized as non-seasonal, five systems showed significant correlation between daily public supply volumes and daily maximum temperature (p<0.001), three showed significant correlation for daily precipitation totals (p<0.05), three showed significant correlation for days since significant precipitation (p<0.05), and three showed significant correlation for the weekday-weekend effect (p<0.01; table 4).

Data Transformations

Prior to fitting a linear regression model to the data, long-term base line trends were estimated and removed from the daily public supply water-use data. A line was fit to the data from 2016 through 2019 for each of the 15 DWSAs to remove any non-stationarity, or multi-year trends, observed in the data (fig. 22). These first 4 years of data were used to “train” the model and the last year of available data (2020) was used to “test” the model. As noted previously, the Washington-Oxford system data have a noticeable step increase between 2019 and 2020. Fitting a linear rate of change to the Washington-Oxford dataset may not have been the most appropriate choice for this system, but it was ultimately used in this study to keep consistency in the analysis across all DWSA systems. The equation used to remove any linear, multi-year trends was of the form

Once the linear multi-year trends were removed from the daily public supply data, daily deviations from mean monthly values were calculated based on the detrended data. To obtain the deviations, mean monthly values were first determined for the detrended daily public supply data (fig. 23) and for daily maximum temperature and daily precipitation totals datasets during the model training period of 2016 through 2019. Mean monthly values were then subtracted from each daily public supply volume, daily maximum temperature, and daily precipitation datapoint to obtain each dataset’s daily deviation from its mean monthly value.

Daily deviations were used as model inputs for multiple reasons. First, the primary study objectives were centered around estimating daily public supply water use and identifying factors that drive daily demand. Therefore, daily fluctuations in public supply water use were isolated from seasonal and annual fluctuations, allowing the potential for influential factors to become more apparent. Second, after studying the variable correlations between the public supply data and predictor variable datasets and analyzing model diagnostics via a review of model residual distribution, using daily deviations from the mean as model inputs showed improved linear correlations and closer adherence to linear regression model assumptions. These transformed daily deviations from the mean data were then used to build a regression model of the form

In addition to the predictor variables discussed previously, daily temperature lagged 1, 2, and 3 days and daily precipitation lagged 1, 2, and 3 days were also included in the initial model development based on previous studies’ findings (Wong and others, 2010; Ahmed and others, 2020). Lagged temperature and precipitation data are defined here as data from either 1, 2, or 3 days prior to the modeled daily public supply water-use day of interest. For example, temperature lagged by 1 day would involve correlating daily public supply water-use data on a given day with temperature data from 1 day prior to that given day. All temperature and precipitation variables used in the model were incorporated as daily deviations from monthly averages, to correspond to the model response variable (daily deviations of public supply).

Selection of Predictor Variables

Backwards stepwise regression selection was used to identify which predictor variables were influential and should be retained in the regression equation for each of the 15 NJAW DWSA systems with daily public supply data. This model selection process involves starting with a regression model that includes all possible predictor variables, and then iteratively removes one predictor variable at a time, until the best fit model is found. The best fit model was determined based on the Akaike information criterion, a metric that is used to estimate the quality of a regression model (Akaike, 1974; R Core Team, 2019; Helsel and others, 2020). The regression variables that were identified as influential from the backwards stepwise regression selection process were then compared amongst the 15 systems, and between the seasonal and non-seasonal groups. The variables that were retained and considered as influential in at least half of the systems in a group (for example, in at least five of the nine systems in the seasonal group, and in at least three of the six systems in the non-seasonal group) were selected for the two regression model forms.

Among the seasonal group, the variables included in the model form were daily maximum temperature, daily maximum temperature lagged by 1 day and 2 days, precipitation lagged by 1 day, 2 days, and 3 days, number of days since significant precipitation, and the weekday-weekend effect. The linear regression equation for the seasonal group is written as

For the non-seasonal group, the variables included in the model form were daily maximum temperature, number of days since significant precipitation, and the weekday-weekend effect. The linear regression equation for the non-seasonal group is written as

The estimated regression coefficients are shown in table 5 for the seasonal group and in table 6 for the non-seasonal group. The regression coefficients represent the correlation factor between each particular predictor variable and the daily public supply water-use deviations. Positive values indicate a positive correlation between the two variables and negative values indicate a negative, or indirect, correlation. The intercept values are empirically derived regression constants.

Autoregressive Integrated Moving Average (ARIMA) Model

One assumption of linear regression models is that the residuals, or errors, are random (Helsel and others, 2020). The residuals from the seasonal and non-seasonal regression model forms based on equations 3 and 4 contained some autocorrelation—when a variable is correlated with itself at different time steps—and were therefore, non-random. This is often the case with time series data as data from any given day are likely to be correlated with data from the previous days. To address this issue of autocorrelation, an autoregressive integrated moving average (ARIMA) model was fitted using the linear regression model residuals with the form

ARIMA models are comprised of an autoregressive component (AR), an integration component (I), and a moving average component (MA) (Hyndman and Athanasopoulos, 2018). These models require three parameters:

For this application, because the data were previously detrended and any long-term trends (or non-stationarity) were removed, there was no differencing applied and, as a result, the d parameter was not relevant in the ARIMA model. The resulting form of the ARIMA model used here is written as

The parameters p and q were chosen based on examination of the partial autocorrelation function (PACF) plots and the autocorrelation function (ACF) plots, as well as use of the “auto.arima” function in the R package “forecast” (version 8.13; Hyndman and others, 2020) which automates the process of identifying the best fit parameters. The PACF and ACF plots help with estimating the values of p and q by showing the amount of autocorrelation in the regression model residuals over time (Hyndman and Athanasopoulos, 2018). Each of the 15 DWSA systems had some variations in the number of parameters but two general equation forms were ultimately chosen—one for the seasonal group and one for the non-seasonal group, to keep consistent with methods developed for this work. The two chosen equation forms were based on averages of the parameters p and q across each group of either seasonal or non-seasonal DWSA systems. Based on the analysis of the PACF and ACF plots in conjunction with the output from the ‘forecast’ R package, only the AR terms were found to be appropriate, and no MA terms were included for both groups. For the seasonal group, the ARIMA model of the form ARIMA(3,0,0) where p=3, d=0, and q=0, was used:

For the non-seasonal group, the model form ARIMA(6,0,0) where p=6, d=0, and q=0, was used:

The ARIMA model coefficient values and whether they were significant at a 95-percent confidence level are provided for the seasonal group (table 7) and non-seasonal group (table 8). The model forms were kept consistent across all systems within each group and therefore were not uniquely specific to each system. For example, all systems in the seasonal group had the same set of predictor variables and number of coefficients included in the model. However, the magnitudes of the model coefficients were unique to each individual system. Because the model forms were not uniquely specific to each system, 84 percent of all the model coefficients across all 15 systems were found to be significant.

There were less terms, defined as the product of the model coefficient and predictor variable, retained in the linear regression models for the non-seasonal systems (3 terms) than for the seasonal systems (8 terms), and more terms retained in the ARIMA models for the non-seasonal systems (6 terms) than for the seasonal systems (3 terms). This indicates that climatic (weather) and weekly variables, generally, are better correlated with daily water use among seasonal systems than non-seasonal systems. Additionally, daily water-use volumes from days prior to the present are better correlated with daily water-use volumes at present, among the non-seasonal systems. In summary, based on the predictor variables examined in this study, daily public supply water use is better predicted by previous days’ water use more so than external variables, such as weather-related factors, for the non-seasonal systems.

Evaluating Regression Models

The primary measures of model error and accuracy used for this study were the root mean squared error (RMSE) and the adjusted coefficient of determination ($R_{adj}^{\text{2}}$). The RMSE can be described by

The RMSE is an absolute measure of model error, based on the difference between the observed and predicted values (Hyndman and Athanasopoulos, 2018). The $R_{adj}^{\text{2}}$ can be described by

The $R_{adj}^2$ is different from the coefficient of determination (R²) in that it also accounts for the number of predictor variables included in the model (Helsel and others, 2020). Adding predictors to a model always increases the R² value because there is an assumption that all predictors explain some of the dependent variable variance. The $R_{adj}^2$ controls against this automatic increase and only increases if new predictor variables improve the model more than would be expected by chance. The $R_{adj}^2$ provides a relative measure of model error and typically ranges from 0 to 1, with 0 indicating 0 percent of the variance in the observed data was explained by the model and 1 indicating 100 percent of the variance was explained by the model. Using a combination of RMSE and $R_{adj}^2$ can provide a better evaluation of the models’ performances (Helsel and others, 2020).

The $R_{adj}^2$ values of the models for each of the 15 DWSA systems are provided in table 9, where the values for the multi-year trends combined with the mean monthly values (column A); with the mean monthly values and modeled daily deviations (column B); and with the mean monthly values, modeled daily deviations, and the ARIMA models (column C) are all presented. The version of the model that included the combination of multi-year trends, mean monthly values, modeled daily deviations, and ARIMA models together showed improved performance and the highest $R_{adj}^2$ values (table 9, column C). Subsequently, results and discussion presented hereafter use this version of the model that is referred to as the linear regression model with autoregressive errors, or LRA model.

In addition to calculating $R_{adj}^2$ values, the models were also visually assessed by plotting and comparing the observed values versus the LRA model values, which were comprised of the multi-year trends, mean monthly values, modeled daily deviations, and ARIMA models (figs. 24 and 25). Where the combination model does well, the points are more tightly aligned with the one-to-one (1:1) reference line, whereas the points are more spread out and scattered for DWSAs for which the model does not perform as well, showing distance from the 1:1 reference line. Generally, the seasonal systems fit closer to the 1:1 reference line than the non-seasonal systems (fig. 25). For some systems, the model underestimates very high values (figs. 24A, 24D, 24G, 24H, and 24I) and overestimates very low values (figs. 24A, 24D, 24G, and 24H), because there are generally fewer extreme high and low data points with which to train the model. A higher $R_{adj}^2$ value indicates less error in the model.

Testing Regression Models on 2020 Data

With regards to using the 2020 data to test the models, because 2020 was the first year of the coronavirus disease 2019 (COVID-19) pandemic, one might expect to see different trends and patterns in water usage. Upon visual inspection, it was not obvious that data from 2020 were different than from previous years, either for daily data from NJAW or monthly data from NJWaTr. It is conceivable that the location and amount of water use changed with the release of certain recommendations from the health department. Stay-at-home guidelines could have shifted water use away from places of employment or recreation because people were advised to stay in their homes and, therefore, their home DWSA. Sanitation guidelines, such as increased hand washing, also potentially increased water use (Kim and others, 2021). For this study, there is insufficient information to determine if the COVID-19 pandemic impacted water usage and more years’ worth of data, in the years preceding and succeeding 2020, would be needed to confirm this theory.

The model performance for each of the 15 case-study systems is shown for the seasonal group (fig. 26) and non-seasonal group (fig. 27). Generally, the models for the seasonal DWSA systems performed better than those for the non-seasonal DWSA systems based on the $R_{adj}^2$ values. The average $R_{adj}^2$ value for the model test year among the seasonal group was 0.78; the average $R_{adj}^2$ value for the same period among the non-seasonal group was 0.25. A possible explanation for why the seasonal DWSA system models perform better might be that the seasonal signal is fairly consistent and predictable from year to year, and the predictor variables used to build the model are mostly weather- or climate-related and thus better correlated with seasonality, which ultimately reduces unexplained variability. It is also possible that there might be additional variables which were not considered in this study that may help to better explain the day-to-day variability in the non-seasonal DWSA systems. However, it may also be that the day-to-day variability is simply random to some degree, and the random variability makes up a larger portion of the data in non-seasonally driven systems because of the lack of a dominant seasonal signal.

Visual comparison of the model training period and test year shows less day-to-day variability in the model estimates during the test year; this is particularly true of the non-seasonal systems during that time (figs. 27D and 27E). The difference in variability between the training period and the test year is due to the “forecasting” or “prediction” method for the ARIMA portion of the models. To estimate daily water use for the 2020 test year, multi-step (or dynamic) ARIMA forecasts were generated. For this estimation method, real or observed data were not used to adjust or correct the regression model estimates as it was assumed the time period of interest is in the future and observed data do not yet exist. In this case, estimating, or “predicting,” water use for 2020 was treated as if it was the future and therefore, no observed data were available to compare to and subsequently correct the model. In contrast to multi-step forecasting of the ARIMA portion of the models, the model predictions from the training period (2016–19) did incorporate the observed daily water-use volumes to correct and adjust the estimates for each day throughout the training period. As an example, when the model was used to estimate daily water use for June 10, 2018, data from 1 through 6 days before (June 4–9, 2018) were incorporated to compute the estimate. This is different from when the model was used to estimate daily water use for a given day in 2020 (test year). For June 10, 2020, for example, December 31, 2019, was the most recent observed data available to adjust the model—almost 6 months prior. The latter method of obtaining model estimates is considered forecasting because it assumes there are no available data for the time period of interest with which to compare and correct the model. With no new data to consider, the forecasted values from the ARIMA portion of the model will converge to zero because of the ARIMA model assumptions (Hyndman and Athanasopoulos, 2018). Ultimately, the predicted variability in daily water use for the 2020 test year is attributed only to external predictor variables (temperature, precipitation, day of the week), and not to autocorrelation, or the portion of variability explained by the ARIMA model, as is by design of the model.

The regression model for the Logan system performs well on the yearly and daily scale (fig. 26H). This system has a noticeable weekly signal where daily public supply water use drops over the weekend and is higher during the weekdays. Because of this regular and predictable pattern, along with an overall seasonal pattern, the regression model is able to capture and predict the day-to-day variations well based upon visual examination and mathematical calculation, with an R² value of 0.76 for the test year (2020). In contrast, the regression model for the Belvidere system does not accurately predict daily public supply water use for 2020 (fig. 27E) based upon visual examination and mathematical calculation, with an R² value of 0.31 for the test year. This result is likely because the Belvidere system does not have a prominent seasonal signal and the day-to-day variability in this system does not appear to be regular or predictable and is seemingly random based on the weather and climatic variables considered in this study. The results from the Washington-Oxford system appear acceptable with an R² of 0.96 (fig. 27F), however it is important to note the exceptional situation for this system. As mentioned previously, there is a large step increase in daily water-use volumes between 2019 and 2020 for this system and explanations for this step increase are unknown (see “Daily Public Supply Water-Use Data” section). Although the same analysis as all other systems was applied here, there is less confidence in the accuracy and meaning in the results because of the large shift in the dataset for the Washington-Oxford system.

The RMSE was also calculated to help evaluate the performance of the models between the training period and the test year (table 10). Because RMSE is an absolute measure of error and the magnitudes of water use vary from system to system, RMSE is not as robust as $R_{adj}^{\text{2}}$ for comparing performance between models (or DWSA systems). Instead, it can provide some insight toward the model error when comparing different time periods from the same model or within each DWSA system. In this case, the training period and the test year were compared.

Anomalous Data Points

As noted previously, in the section titled “Daily Public Supply Water-Use Data,” there were a few observed data points that were outside the typical range of values for a given DWSA system. The same regression analysis was performed on the datasets with these data points removed to evaluate the effect of these anomalous data on the regression models (fig. 28). There were five systems that had visually identifiable anomalous data points (Burlington, Strathmere, and Bridgeport in figure 26; Camden and Frenchtown in figure 27). When the regression analysis was applied with and without these anomalous data points, the differences in $R_{adj}^2$ values ranged from 0.001 to 0.15 with a mean difference of 0.04 (table 11). The differences in RMSE ranged from 0.001 to 0.05 Mgal/d, with a mean difference of 0.01 Mgal/d (table 12). These findings indicate that including the few anomalous data points does not drastically affect the models and their performance (tables 11 and 12). Thus, retaining anomalous data points in any additional datasets that may be used to test the models may not have a sizable effect on the results.

Table 11.    

The adjusted coefficient of determination ($R_{adj}^{\text{2}}$) values for the New Jersey American Water drinking water service area (DWSA) systems that had anomalous data, shown for the training period (2016–19) and the test year (2020) with the anomalous data removed and retained for comparison.

[PWSID, public water system identification number]

Table 11.    The adjusted coefficient of determination ( R a d j 2 ) values for the New Jersey American Water drinking water service area (DWSA) systems that had anomalous data, shown for the training period (2016–19) and the test year (2020) with the anomalous data removed and retained for comparison.

PWSID	DWSA system	Training period		Test year
		$R_{adj}^{\text{2}}$ without anomalous data	$R_{adj}^{\text{2}}$ with anomalous data	$R_{adj}^{\text{2}}$ without anomalous data	$R_{adj}^{\text{2}}$ with anomalous data
NJ0327001	Burlington1	0.78	0.78	0.73	0.71
NJ0327001	Camden1	0.42	0.34	0.11	0.08
NJ0511001	Strathmere	0.81	0.78	0.87	0.87
NJ0809001	Bridgeport	0.44	0.42	0.58	0.43
NJ1011001	Frenchtown	0.49	0.49	0.01	0.00

¹

System represents a portion of the Western Division system and therefore its data represent a portion of the Western Division system's water-use data available in NJWaTr.

Table 12.    

The root mean squared errors for the New Jersey American Water drinking water service area (DWSA) systems that had anomalous data, shown for the training (2016–19) and test (2020) periods with the anomalous data removed and retained for comparison.

[PWSID, public water system identification number; RMSE, root mean squared error]

Table 12.    The root mean squared errors for the New Jersey American Water drinking water service area (DWSA) systems that had anomalous data, shown for the training (2016–19) and test (2020) periods with the anomalous data removed and retained for comparison.

PWSID	DWSA system	Training period		Test year
		RMSE without anomalous data	RMSE with anomalous data	RMSE without anomalous data	RMSE with anomalous data
NJ0327001	Burlington1	0.520	0.520	0.834	0.868
NJ0327001	Camden1	0.106	0.122	0.181	0.227
NJ0511001	Strathmere	0.020	0.022	0.023	0.023
NJ0809001	Bridgeport	0.017	0.018	0.021	0.025
NJ1011001	Frenchtown	0.015	0.017	0.044	0.054

¹

System represents a portion of the Western Division system and therefore its data represent a portion of the Western Division system's water-use data available in NJWaTr.

Currently, the models do not estimate water use on the days with anomalous data well and may not elucidate anomalous data points found in additional datasets (for other time periods, beyond the training period and the test year, or for other DWSA systems). The reasons behind some of these data points, such as water main breaks (NJAW, written commun., 2021), are usually unpredictable. The goal of these regression models was to estimate typical daily public supply water use, so it is expected that the models may not accurately estimate these atypical days and events. Additionally, not all predictor variables and information that might be relevant in predicting these events, such as the age of the water supply system and its infrastructure, were included in this study. Even though the models do not predict these anomalous data points well, the model performance statistics are mostly unchanged even with these points removed owing to the small number of anomalous data points relative to the entire dataset. If a dataset were to have many more data points with uncharacteristic volumes of water use, the anomalous data points could have a larger impact on the overall model and its performance.

Accuracy and Limitations of Regression Equations

The primary measures of model error and accuracy used in this study were the RMSE and $R_{adj}^2$. Comparing each statistic between the training period and the test year, the models’ performance tends to drop slightly during the test year. The average increase in RMSE for the test year is 0.13 Mgal/d. The average decrease in $R_{adj}^2$ for the test year is 0.18 Mgal/d. This is to be expected as the test data were treated as “new” data for the model and were not involved in training or developing the model. Instead, the model had to “predict” or extrapolate the data from the test year.

These models were created for the 15 DWSA systems in New Jersey for which daily public supply water-use data were provided and are therefore only appropriate for estimating water use within those 15 specific systems. There are likely commonalities in factors affecting daily public supply water use across DWSA systems in New Jersey. However, more daily data (from additional systems and additional years) are needed to better identify those commonalities and the extent to which these equations developed here are applicable to other systems. Additionally, there are likely influential factors beyond those studied here that could potentially improve the accuracy of daily public supply water-use estimates. The factors included in this study are commonly used factors in estimating water use, however there may be additional factors specific to New Jersey that could be suitable (Zhou and others, 2000; Eslamian and others, 2016; Opalinski and others, 2019; Ahmed and others, 2020). Finally, having more years’ worth of daily water-use data to train the models would likely improve the model performance, particularly for identifying any long-term trends that may be present in the data.

Developing and Using the Generalized Regression Models

This method for monthly-to-daily disaggregation would likely be more useful if it could be applied to systems beyond the 15 DWSA systems for which daily water-use data were acquired and for which the daily regression models were developed. To increase the applicability of this tested disaggregation method, a single averaged regression equation was developed to represent all seasonal systems, and a second averaged regression equation was developed to represent all non-seasonal systems. These averaged equations were developed based on the 15 case-study systems’ equations, but the intent was for them to be applicable for any system in the state, which assumes the 15 systems were representative of all DWSA systems in the state. However, having more daily data from more systems to further develop the averaged equations would likely improve their accuracy.

To compute the averaged regression equations, also referred to as the generalized regression equations, the daily public supply water-use data were first normalized by the populations served for each system. Then, the mean monthly values were calculated and subtracted from the data to obtain daily deviations. The normalized daily deviations public supply data were then averaged across each system grouping (seasonal and non-seasonal). For example, the normalized, daily deviations public supply data were averaged across all nine seasonal systems for each day in the period of interest (2016–20) resulting in one set of averaged (and normalized) daily deviations public supply data for systems demonstrating seasonal water-use patterns. This averaging was done for the climatic variables as well. For example, the daily maximum temperatures were averaged at each time step across all nine seasonal systems to obtain one averaged, daily maximum temperature.

Next, the same regression development method was applied as outlined above. Backward stepwise regression selection was used to select which predictor variables were best suited to predict the averaged, normalized daily public supply water-use deviations (eqs. 11 and 12). Once the best suited variables were identified, the regression model coefficients were estimated using R statistical software (R Core Team, 2019). This resulted in two sets of averaged (and normalized) model coefficients, one for seasonal systems and one for non-seasonal systems (table 13).

ARIMA models for the generalized seasonal and non-seasonal groups were also developed using the same method from the individual daily ARIMA models (table 14). The seasonal system equation is of the form:

The non-seasonal system equation is of the form:

To test the accuracy of the generalized equations on a system for which there were monthly NJWaTr data but no daily data, the monthly data were first normalized by that system’s served population to scale the water-use volumes to the same magnitude.

Next, the generalized regression equations were tested to determine if combining these equations with disaggregated monthly-to-daily NJWaTr data could produce reasonable daily water-use estimates for DWSA systems beyond the 15 highlighted in this study. The generalized equations were combined with the disaggregated monthly NJWaTr data using the same process discussed at the beginning of this section, but instead of using system-specific regression equations, the two generalized models (seasonal and non-seasonal) were used. The final step was to multiply the water-use volumes and their respective DWSA system populations served to rescale the data to the system’s original magnitudes. The estimated daily public supply water-use volumes were then compared to the NJAW daily public supply water-use data from 2016 through 2020 for the 15 systems with daily data (figs. 32 and 33). Overall, the estimates from the generalized models were similar to the system-specific models, with an overall average 6-percent increase in RMSE, which indicates this method is a reasonable way to estimate daily public supply given monthly data and accompanying climate information. In general, as observed previously with the system-specific regression equations the generalized seasonal equations performed better than the generalized non-seasonal equations and had higher adjusted R² values (figs. 32 and 33), especially when the daily data align well with the monthly data from NJWaTr. There are limitations in using these generalized regression equations for systems where only partial data were provided, and therefore the data did not represent the whole system. Limitations also exist for systems where the aggregated daily data did not match the monthly totals in NJWaTr. These limitations are described in more detail in the “Limitations of Generalized Regression Models” section.

Testing Generalized Regression Models on New DWSA Systems

As a final test of these generalized regression equations, three seasonal DWSA systems for which daily data were not available were selected from the NJWaTr database. The three systems were chosen based on their similarity to three seasonal NJAW systems. The first system was Long Beach–North Beach (NJ1517003; table 15, Comparison Group A), which was chosen because of its similarity to the Strathmere system: they both have smaller population-served (Q2) values, higher income and property values, non-urban residential density, and coastal geography. The second system chosen was Margate City (NJ0116001; table 15, Comparison Group B), which was similar to the Ocean City system: they both have larger population-served (Q4) values, higher property values, urban residential density, and coastal geography. The third system chosen was Ringwood (NJ1611002; table 15, Comparison Group C), because of its similarity to the Harrison system: they both have larger population-served (Q3) values, higher income and property values, non-urban residential density, and non-coastal geography.

Because the seasonal generalized regression models perform better than the non-seasonal generalized regression models discussed in this study, the focus was to test the generalized seasonal regression model, as there are likely additional influential factors to be explored with the non-seasonal model and modeled results of the non-seasonal systems showed more errors and less correlation. With these three test DWSA systems, the monthly NJWaTr data were disaggregated to an average daily total using the simple disaggregation method of dividing by number of days in the month and the generalized regression estimates (deviations from the mean monthly values) were added to include daily fluctuations. These final estimates were then compared to the disaggregated NJWaTr data, as there were no daily data to compare to for these three systems (fig. 34). The RMSEs were used to assess the modeled estimates in comparison to the disaggregated NJWaTr data and ranged from 0.006 to 0.17 Mgal/d (fig. 34). It is difficult to compare monthly and daily datasets; however, the primary purpose of this exercise was to determine if it is possible to obtain reasonable daily water-use volumes in the absence of daily data. The final estimates may be biased (high or low) because they are being assessed against data that were used to generate the estimates.

Limitations of Generalized Regression Models

In addition to the limitations and extent of accuracy discussed in the section titled “Accuracy and Limitations of Regression Equations” with the 15 system-specific regression equations, there are additional constraints to consider that are introduced when using the generalized regression models and incorporating the NJWaTr dataset of monthly values. Combining the monthly NJWaTr datasets into the water-use estimates creates some additional uncertainty because the water-use data used in the estimations comes from two separate sources: one from NJAW’s daily dataset and one from NJWaTr’s monthly dataset. Both datasets may have some uncertainty and error that are likely not the same across, and potentially within, each dataset. Because these datasets may have different sources of uncertainty and types of errors, there is potential for the summed daily values not to match the monthly values, as was found in several systems. Having a better understanding of those systems whose daily values do not sum to monthly values could be valuable to improve model results and outcomes.

For the NJAW dataset, the level of quality assurance and quality control is unknown; for example, the specifics of meter calibration or other types of verification are unknown. Some examples of uncertainty and factors affecting accuracy were presented in the discussion on anomalous data points (in the section titled “Anomalous Data Points”) that were observed in several systems with daily data from NJAW. Data from NJWaTr were checked using a quality assurance and quality control procedure in which monthly data are checked against permitted allocation limits and checked against historical data for each site for its period of record (Shourds, 2020). Another source of uncertainty or errors inherent to the NJWaTr dataset may result from the fact that these public supply data come from many different purveyors. Furthermore, each system within a public supply purveyor could have its own individual sources of error or conditions specific to that system that may affect the quality and accuracy of the data.

Another source of uncertainty arises from averaging regression equation coefficients. Each model was initially developed for a specific DWSA system and therefore had unique coefficient values based on the correlation between a particular system’s daily water use and various predictor variables (like temperature and precipitation). By averaging the coefficients, any interpretation on the correlations between variables is not recommended. Rather, the primary purpose of averaging the coefficients was to determine if reasonable estimates of water use could be made on systems beyond the 15 systems for which daily data were acquired. This method of averaging regression coefficients and combining the outputs with NJWaTr data was tested for this study without any preconceived knowledge of what the results would be.

More data and further analysis would be needed to determine if this is a viable and accurate method for estimating daily public supply water use for any DWSA system in New Jersey. More daily data would allow for further model output verification. Additionally, exploring other factors or variables that may influence daily public supply water use, beyond those that were considered in this study, could also improve model outcomes.