Abstract

Wastewater flow forecasting is key for proper management of wastewater treatment plants (WWTPs). However, to predict the amount of incoming wastewater in WWTPs, wastewater engineers face challenges arising from numerous complexities and uncertainties, such as the nonlinear precipitation-runoff relationships in combined sewer systems, unpredictability due to aging infrastructure, and frequently inconsistent data quality. To address such challenges, a time series analysis model (i.e., the autoregressive integrated moving average, ARIMA) and an artificial neural network model (i.e., the multilayer perceptron neural network, MLPNN) were developed for predicting wastewater inflow. A case study of the Barrie Wastewater Treatment Facility in Barrie, Canada, was carried out to demonstrate the performance of the proposed models. Fifteen-minute flow data over a period of 1 year were collected, and the resampled daily flow data were used to train and validate the developed models. The model performances were examined using root mean square error, mean absolute percentage error, coefficient of determination, and Nash–Sutcliffe efficiency. The results indicate that both models provided reliable forecasts, while ARIMA showed a slightly better performance than MLPNN in this case study. The proposed models can provide useful decision support for the optimization and management of WWTPs.

INTRODUCTION

Influent flow rate is an important parameter for the operation of wastewater treatment plants (WWTPs). Not only could it provide guidance for wastewater treatment operation to achieve stable effluent quality, it could also help identify capacity risks and develop alternatives to upgrade existing infrastructure. WWTP managers and operators are in urgent need of efficient practical tools to predict incoming wastewater flow, in order to optimize the treatment processes and thus reduce environmental risks and improve treatment efficiency.

Previously, a number of wastewater forecasting models have been developed based on the simulation of wastewater collection systems (Ashley et al. 1999; Benedetti et al. 2013). For example, Schmitt et al. (2004) used a drainage simulation model based on hydraulic flow routing to develop a management tool for urban drainage systems. Lindqvist et al. (2005) used a physical model and an autoregressive method to predict influent flow to the Rya WWTP in Sweden. However, such models have several limitations. For example, the physical processes that affect wastewater inflow, such as rainfall, snowmelt, infiltration and runoff, are typically too complex to simulate. In addition, wastewater flow is affected by domestic water use, which is difficult to be integrated into urban drainage models. Moreover, the aging of wastewater infrastructure (such as pipes and sewer connections) makes it even more difficult to develop physical models to predict wastewater inflow, due to the changes in facilities (Tan et al. 1991; Ashley et al. 1999).

To address the challenges in the simulation of such complex systems, numerous data-driven techniques have been developed. Due to their effectiveness and efficiency, approaches such as the autoregressive (AR) model, artificial neural network (ANN), genetic algorithm, and multivariate analysis have been widely applied. Particularly, several studies of inflow forecasting at WWTPs using data-driven techniques have been reported (Djebbar & Kadota 1998; El-Din & Smith 2002; Wei et al. 2013). For example, a time series model was developed by Kim et al. (2006) to forecast flow rate and chemical components in the wastewater influent, which yielded good results for 1- to 7-day forecasting. Fernandez et al. (2009) developed a fuzzy neural network model and obtained forecasting results with an average error below 10%. Kim et al. (2016) used k-nearest neighbor method to predict the influent flow rate and four water qualities, and finally achieved good predicted results. These studies provided proof-of-concept for the application of data-driven techniques in wastewater forecasting (Campisano et al. 2013). However, no successful practical application of data-driven techniques for wastewater forecasting in Canada has been reported to date.

In this study, a time series analysis model (i.e., the autoregressive integrated moving average, ARIMA) and an ANN model (i.e., the multilayer perceptron neural network, MLPNN) will be developed for inflow forecasting. The proposed models will be applied to the Barrie Wastewater Treatment Facility (WwTF) in Barrie, Ontario, Canada. This plant collects sanitary wastewater from the city's residential, commercial, and industrial water users. It also experiences significant inflow and infiltration during rainfall events as a result of the downspouts and sump pumps illegally connected to the sanitary system. This study entails the following key objectives: (1) analyze the temporal patterns of influent flow at the Barrie WwTF and develop a procedure for wastewater characterization; (2) develop an ARIMA model for predicting daily wastewater inflow at the Barrie WwTF; and, finally, (3) develop an MLPNN model for daily wastewater inflow prediction and compare it with the proposed ARIMA model.

METHODOLOGY

Modeling development

Autoregressive integrated moving average

The ARIMA model was developed by Box and Jenkins in 1976 (Box & Jenkins 1976). Soon after it was developed, ARIMA was widely applied in hydrological modeling (Machiwal & Jha 2012). Generally, the ARIMA model often produces satisfactory results, and it can describe the variation in flow rate on a daily, monthly and annual time step (Al-Asheh et al. 2007). An ARIMA model consists of three components, which are the AR process, the moving average (MA) process, and the integrated component (I). The AR process assumes that each observation is a linear combination of prior observations and a random error component . The AR equation can be written as: 
formula
(1)
where are the AR model parameters. The highest value of p, for which , is the order of the AR process. Thus, the AR process of order p can be denoted by AR(p).
The MA process assumes that each observation is a linear combination of prior error components, and a random error component . The MA equation can be written as: 
formula
(2)
where are the MA model parameters. The highest value of q, for which , is the order of the MA process. Thus, the MA process of order q can be denoted by MA(q).
The I component is for addressing the non-stationarity of the time series. A time series is stationary if its statistical properties (e.g., mean, variance, autocorrelation) are constant over time (Machiwal & Jha 2012). Since the data series in ARIMA needs to be stationary, a non-stationary time series needs to be transformed through differencing. Differencing eliminates the trend and seasonality of time series data by calculating the difference between consecutive observations (Equation (3)). Sometimes the data should be differenced twice to obtain stationarity, which is known as second order differencing (Equation (4)). The number of transformations that are required to obtain stationarity is denoted by d. A process which has dth order differencing transformations is called an integrated process of order d, and it can be denoted as I(d). 
formula
(3)
 
formula
(4)

Thus, an ARIMA model with the order of (p, d, q) is a combination of AR(p), MA(q), and I(d). For example, an ARIMA (p, 0, 0) process with is an AR(p) process, whereas an ARIMA (0, 0, q) process with is an MA(q) process. In ARIMA, it is assumed that the random error component is a randomly drawn series from a normal distribution, which is a white noise process (Contreras et al. 2003).

Multilayer perceptron neural network

ANN is a machine learning technique which could be used to build predictive models by extracting information from past experience. It shows significant improvements of model flexibility compared to traditional predictive models (Novotny et al. 1991). Research on ANN began in 1943 (McCulloch & Pitts 1943) and developed rapidly in many fields in the following decades (Fletcher & Goss 1993; Nielsen et al. 1997; Ordieres et al. 2005).

The term ANN reflects its analogy to the human nervous system, and an ANN consists of several simple computational elements called ‘artificial neurons’. Because of its simpler topology compared to many other ANN models, the multilayer feedforward neural network (MLFFNN) is widely used in science and engineering to solve the problems involving classification, regression, and time series analysis (Rafiq et al. 2001). In MLFFNNs, neurons are sorted in layers, where each layer consists of one or more neurons: an input layer, an output layer, and n hidden layer(s), where n could be adjusted according to the model performance. Each neuron in one layer has only direct connection to a neuron in the next layer by a connected weight, which reflects the strength of the relationship between two neurons. The weights are initialized to small random numbers and updated to minimize the prediction error through an optimization process.

The ANN is fed by a series of input-output pairs and is trained to reproduce the outputs. The learning process of a feedforward network occurs in the perceptron, which consists of a single neuron with adjustable weights and an activation function. When an MLFFNN has perceptrons with two or more trainable weight layers, it might also be called an MLPNN. In this study, a back-propagation algorithm is used to train the MLPNN. This optimization algorithm minimizes the error between observations and predictive outputs by adjusting neuron weights. Essentially, it evaluates the gradient of an error function, e.g., mean square error, to adjust the weights for a feedforward network (Jain et al. 1996).

In this study, an MLPNN model was developed for wastewater flow forecasting. The inputs were prior daily flow rate data, daily rainfall, daily snow, daily snow on the ground, 1–7 days prior precipitation accumulations, the index of day of the week, and daily air temperature (max, min, and mean). Seven binary variables were used to indicate the day of the week. For example, if it is a Monday, the corresponding variable takes a value of 1 and all other six variables are assigned a value of 0.

A three-layer MLPNN was selected as the best model structure, based on a series of exploratory experiments. The exploratory experiments manually searched and tested different numbers of hidden layers, neurons, and epochs. The MLPNN used in this study consists of one input layer with 21 input variables, one output layer with one output variable, and one hidden layer with 20 neurons (Figure 1). The epoch used in this study is 100.

Figure 1

Architecture of the developed MLPNN model.

Figure 1

Architecture of the developed MLPNN model.

Evaluation of model performance

In order to examine the model accuracy, different evaluation criteria, such as root mean square error (RMSE), mean absolute percentage error (MAPE), coefficient of determination (R2), and Nash–Sutcliffe efficiency (NSE), were used in this study.

RMSE and MAPE are two statistical indexes which are widely used for measuring the errors between predictions and observations. A model with low RMSE and MAPE is preferred. R2 with a range of [0, 1] is a statistical measure of how close the predictions are to the real observations data (Bennett et al. 2013). An R2 value close to 1 indicates the fitness of the predictions to the real observed data. NSE, defined by Nash & Sutcliffe (1970), has a range of . If NSE < 0, the model performs worse than its mean value; an NSE value of 1 indicates a perfect fit between the observations and the predictions.

Study area

Barrie wastewater treatment facility

The Barrie WwTF in Barrie, Ontario, Canada, is the largest WwTF on the shores of Lake Simcoe (Figure 2). According to a 2016 census (Statistics Canada 2017), the City of Barrie has a total area of 898.02 km2 and a population of 141,434, which made the city the 34th largest city in Canada in terms of population. The existing sanitary infrastructure system in the city includes 530 km of sewers, 9 km of force mains, 15 active pumping stations (during the duration of the study), and one sewage treatment plant (i.e., the Barrie WwTF).

Figure 2

Locations of the Barrie WwTF and meteorological stations (Basemap: ESRI Light Grey Canvas. Source: Esri, HERE, DeLorme, MapmyIndia, © OpenStreetMap contributors, and the GIS user community (Esri 2011)).

Figure 2

Locations of the Barrie WwTF and meteorological stations (Basemap: ESRI Light Grey Canvas. Source: Esri, HERE, DeLorme, MapmyIndia, © OpenStreetMap contributors, and the GIS user community (Esri 2011)).

The storm sewer system and sanitary sewer system are separate in the City of Barrie. The Barrie WwTF is designed to collect sanitary wastewater from the city's residential, commercial, and industrial water users. However, because of the downspouts and sump pumps illegally connected to the sanitary system, it does experience significant inflow and infiltration during rainfall events. During some high flow events, the City of Barrie surcharges the collection system by ‘peak shaving’, which would cap the amount being pumped into the plant and allow excess flow to remain in the sewer until after the peak flow has passed. The influent flow being reported in the paper is the amount being pumped from the wet well to the treatment system. The design capacity of the facility is 76,000 m3/day and the design peak flow is 156,000 m3/day. The Barrie WwTF provides a tertiary treatment system to ensure that the effluent quality meets the water quality standards of Lake Simcoe before sending it back into the lake.

Data collection

The influent flow rates were measured on a 15-minute basis from November 1, 2015 to October 31, 2016. There are four missing flow data points on March 13 from 2:00 to 3:00 am due to the daylight savings time; the total number of flow rate data obtained is 35,132.

Daily meteorological data in the same time period, including maximum, minimum and mean air temperature, total rainfall, total snow, total precipitation, and snow on the ground were obtained from Environment Canada's National Climate Data and Information Archive. Among the five meteorological stations near the WwTF, two stations (Barrie Landfill 44.39°N, 79.74°W, and Shanty Bay 44.40°N, 79.63°W) were chosen for collecting meteorological data (Figure 2). These two stations were selected because they are located closer to Barrie WwTF and had more meteorological data available than the other three stations. To address the spatial variability of meteorological conditions, the average value of data from these two stations was used to describe the meteorological conditions at the Barrie WwTF.

Data pre-processing

In order to improve the data quality, outliers should be eliminated (Verma et al. 2013; Wei et al. 2013). A robust outlier detection technique known as Hampel X84 was used in this study. A simple version of Hampel X84 considers any data point that is more than 1.4826x median absolute deviation (MAD) away from the median as outliers. The MAD can be calculated by Equation (5), and x is the number of standard deviations away from the mean (Hellerstein 2008). For outlier detection within this study, the values outside the range of C ± 4.4478 MADs (i.e., x = 3) were considered as outliers. 
formula
(5)
where ki are samples, , n is the total number of samples, and C = median {ki}.

Considering the flow rate variations in different seasons (e.g., ice melting in the spring and storm events in the summer), the 1-year data were partitioned to four time periods (four seasons), and outliers were identified on a season-by-season basis. Table 1 provides the upper- and lower-bound values, number, and proportion of outliers in each season. The percentages of outliers are 0.59%, 0.44%, 0.08%, and 0.00% in the spring, summer, fall, and winter, respectively. The outliers accounted for 0.28% of the total raw data and were removed before analyzing the flow pattern and developing the predictive models. The removed outliers and missing meteorological data were filled in by linear interpolation (Junninen et al. 2004).

Table 1

Range, amount and proportion of outliers

 C + 4.4478MADs (m3/hr) C − 4.4478MADs (m3/hr) Total amount Outlier amount (>) Outlier ratio (>) Outlier amount (<) Outlier ratio (<) Outlier ratio (total) 
Spring 4,221.07 1,161.99 8,828 50 0.57% 0.02% 0.59% 
Summer 3,528.44 1,258.62 8,832 30 0.34% 0.10% 0.44% 
Fall 3,733.94 835.06 8,736 0.08% 0.00% 0.08% 
Winter 3,843.89 1,073.88 8,736 0.00% 0.00% 0.00% 
Year – – 35,132 87 0.25% 11 0.03% 0.28% 
 C + 4.4478MADs (m3/hr) C − 4.4478MADs (m3/hr) Total amount Outlier amount (>) Outlier ratio (>) Outlier amount (<) Outlier ratio (<) Outlier ratio (total) 
Spring 4,221.07 1,161.99 8,828 50 0.57% 0.02% 0.59% 
Summer 3,528.44 1,258.62 8,832 30 0.34% 0.10% 0.44% 
Fall 3,733.94 835.06 8,736 0.08% 0.00% 0.08% 
Winter 3,843.89 1,073.88 8,736 0.00% 0.00% 0.00% 
Year – – 35,132 87 0.25% 11 0.03% 0.28% 

RESULTS AND DISCUSSION

Characteristics of the wastewater inflow

The 15-minute wastewater inflow data at Barrie WwTF were collected from November 1, 2015 to October 31, 2016. The flow rate varied from 1,056 m3/hr to 4,206 m3/hr, with an average rate of 2,335 m3/hr (outliers excluded). The diurnal patterns in each season and each month, as well as each day of the week are presented in Figure 3. Figure 3(a) indicates similar variation patterns among four seasons. In all four seasons, the flow first declined gradually from 10:00 pm and dropped to the minimum at around 6:00 am. The flow then increased rapidly until 11:00 am, and then decreased slightly at 6:00 pm. After 6:00 pm, the flow increased to a second peak at 10:00 pm. The highest inflow recorded was 4,206 m3/hr in the spring due to snow melting, whereas the lowest inflow recorded was 1,056 m3/hr in the fall due to the dry weather. The hourly flow range in each season (Figure 4) shows that the uncertainty of flow in spring is the highest among the four seasons, which could be a result of the complexity of snow melting processes. The characteristics of flow ranges are similar in the fall and winter, which implies a relatively stable inflow during these seasons. In the summer, the flow ranges are mostly narrow, with only a few variations, which might be attributed to summer storm events.

Figure 3

Diurnal patterns of flow rate in each season of the year (a), each month of the year (b), each day of the week (c).

Figure 3

Diurnal patterns of flow rate in each season of the year (a), each month of the year (b), each day of the week (c).

Figure 4

Diurnal flow rate change in spring (a), summer (b), fall (c), and winter (d).

Figure 4

Diurnal flow rate change in spring (a), summer (b), fall (c), and winter (d).

The monthly diurnal changes of flow rate are shown in Figure 3(b). Similar to the seasonal patterns shown in Figure 3(a), the highest flow levels occurred in March and April whereas the lowest level was observed in October. The hourly flow range (1,460–4,206 m3/hr) in March is much broader than that in April and May (1,297–3,940 and 1,180–3,565 m3/hr, respectively), which indicates that the large fluctuations during spring mainly took place in March. Based on the synchronous meteorological data of total precipitation, total snow, and air temperature, the high and unstable flow in March resulted from the influx of melting snow when the air temperature was above zero. Compared to March, the flow range in October is narrow and stable (1,056–3,697 m3/hr), which indicates a relatively stable flow during the fall.

Moreover, the diurnal pattern of flow rate in each day of the week was analyzed and shown in Figure 3(c). The characteristics of the influent flow rate in weekdays are similar. The influent flow rates during weekdays decreased gradually from 10:00 pm to 6:00 am and dropped to a minimum of 1,540 m3/hr. Then the flow rates increased rapidly from 6:00 am to around 10:00 am, and stayed at a relatively stable level with a slight decline until 6:00 pm. After 6:00 pm, they increased again until around 10:00 pm and reached the maximum of 2,808 m3/hr, except for Friday with a maximum of 2,739 m3/hr. Compared to the flow characteristics of weekdays, the inflow rates presented a different pattern during weekends. There is an obvious time lag of inflow rate between weekends and weekdays (Figure 3(c)). The flow rates during the weekend decreased gradually from midnight to 7:00 am and dropped to a minimum of 1,503 m3/hr, then began to increase to reach a maximum of 3,116 m3/hr around 12:00 noon. After 12:00 noon, on Saturday the flow rate decreased until the end of the day, while on Sunday flow rate continued to drop until 6:00 pm, then increased until 10:00 pm. As we mentioned before, the Barrie WwTF collects sanitary wastewater from the city's residential, commercial, and industrial water users. Generally, the residential water consumption peaks at 8 am on weekdays and between 10 am and noon on weekends (Buchberger & Wells 1996), which shows the same pattern as the wastewater inflow rate. Therefore, it is indicated that the time lag of flow rates and different peak values between weekdays and weekend were caused by the changes in residential water consumption.

Results of the ARIMA model

Model configuration and parameter estimation

In order to check the stationarity of the raw data, the autocorrelation function (ACF) and partial autocorrelation function (PACF) (Little 2013) of the data were calculated, which showed that the flow rate was not stationary. Since the ACF and PACF of the flow data with a 15-minute interval showed that the time series contained too much variation, the time series data were resampled on a daily basis. However, the ACF and PCAF values showed that the daily time series was also not stationary.

Therefore, the first-order differencing was carried out to ensure stationarity of the daily flow data, and the differenced data are presented in Figure 5(a). In order to test the stationarity of the first-order differenced data, the ACF and PACF (Little 2013) of the first-order differenced data were calculated (Figure 5(b) and 5(c)). After the one-step differencing, the ACF and PACF decay quickly and stay within the confidence interval, which indicates that the differenced data is stationary. Thus, the d value was determined to be 1.

Figure 5

Sequence plot (a), ACF (b), and PACF (c) of first-order differenced daily flow.

Figure 5

Sequence plot (a), ACF (b), and PACF (c) of first-order differenced daily flow.

Grid search, which is a process to determine the optimal combination of model hyperparameters, was used to test different combinations of p (0, 1, 2, 3, 4, 5, and 6) and q (0, 1, 2, and 3) with d = 1 to find the ARIMA model with the best fit (Table 2). The dataset was split into two parts: (1) a training set including 66% of the data and (2) a testing set including the remaining 34%. According to Table 2, in the training set, ARIMA (5,1,1) has the lowest RMSE of 109.03 and the highest NSE of 0.776 compared to other tested models. Thus, the ARIMA (5,1,1) was determined to be the most suitable model for inflow rate forecasting in this study.

Table 2

The training and testing results of suggested ARIMA models

ARIMA (p, d, q) Training
 
Testing
 
RMSE MAPE (%) R2 NSE RMSE MAPE (%) R2 NSE 
(0, 1, 1) 114.55 3.050 0.763 0.752 76.05 2.626 0.586 0.563 
(0, 1, 2) 111.20 2.916 0.772 0.767 72.75 2.507 0.613 0.600 
(0, 1, 3) 111.72 2.948 0.771 0.765 72.31 2.493 0.620 0.605 
(1, 1, 0) 115.94 3.042 0.760 0.747 78.37 2.654 0.573 0.536 
(1, 1, 1) 111.62 2.881 0.768 0.765 74.03 2.552 0.595 0.586 
(2, 1, 0) 111.77 3.006 0.773 0.765 73.67 2.559 0.609 0.590 
(2, 1, 1) 111.06 3.011 0.777 0.768 72.95 2.517 0.620 0.598 
(3, 1, 0) 111.86 3.017 0.774 0.764 73.77 2.568 0.611 0.589 
(3, 1, 1) 113.28 3.012 0.768 0.758 72.71 2.513 0.621 0.600 
(4, 1, 0) 110.48 2.961 0.777 0.770 72.28 2.523 0.622 0.605 
(4, 1, 1) 109.05 2.880 0.778 0.776 71.86 2.530 0.619 0.610 
(5, 1, 0) 110.17 2.908 0.777 0.771 72.68 2.564 0.617 0.601 
(5, 1, 1) 109.03 2.879 0.779 0.776 71.64 2.519 0.621 0.612 
(6, 1, 0) 109.58 2.889 0.779 0.774 71.11 2.497 0.630 0.618 
(6, 1, 1) 109.17 2.862 0.780 0.775 71.63 2.500 0.625 0.612 
ARIMA (p, d, q) Training
 
Testing
 
RMSE MAPE (%) R2 NSE RMSE MAPE (%) R2 NSE 
(0, 1, 1) 114.55 3.050 0.763 0.752 76.05 2.626 0.586 0.563 
(0, 1, 2) 111.20 2.916 0.772 0.767 72.75 2.507 0.613 0.600 
(0, 1, 3) 111.72 2.948 0.771 0.765 72.31 2.493 0.620 0.605 
(1, 1, 0) 115.94 3.042 0.760 0.747 78.37 2.654 0.573 0.536 
(1, 1, 1) 111.62 2.881 0.768 0.765 74.03 2.552 0.595 0.586 
(2, 1, 0) 111.77 3.006 0.773 0.765 73.67 2.559 0.609 0.590 
(2, 1, 1) 111.06 3.011 0.777 0.768 72.95 2.517 0.620 0.598 
(3, 1, 0) 111.86 3.017 0.774 0.764 73.77 2.568 0.611 0.589 
(3, 1, 1) 113.28 3.012 0.768 0.758 72.71 2.513 0.621 0.600 
(4, 1, 0) 110.48 2.961 0.777 0.770 72.28 2.523 0.622 0.605 
(4, 1, 1) 109.05 2.880 0.778 0.776 71.86 2.530 0.619 0.610 
(5, 1, 0) 110.17 2.908 0.777 0.771 72.68 2.564 0.617 0.601 
(5, 1, 1) 109.03 2.879 0.779 0.776 71.64 2.519 0.621 0.612 
(6, 1, 0) 109.58 2.889 0.779 0.774 71.11 2.497 0.630 0.618 
(6, 1, 1) 109.17 2.862 0.780 0.775 71.63 2.500 0.625 0.612 

The lowest RMSE, lowest MAPE, highest R2, and highest NSE values are marked bold.

Model validation

As shown in Table 2, the RMSE, MAPE (%), R2, and NSE for the testing set are 71.64, 2.519, 0.621 and 0.612, respectively, which indicates ARIMA (5,1,1) is an acceptable model. A comparison between the observed and predicted time series is presented in Figure 6(a), indicating a good fit between the observations and their predictions. The scatter plot (Figure 7(a)) also indicates satisfactory forecasting results, since most of the points are close to the diagonal line except for some extreme high values.

Figure 6

Comparison of predicted and observed flow rate values: ARIMA model (a) and MLPNN model (b).

Figure 6

Comparison of predicted and observed flow rate values: ARIMA model (a) and MLPNN model (b).

Figure 7

Scatter plot of wastewater flow rate observations versus predictions by ARIMA model (a) and MLPNN model (b).

Figure 7

Scatter plot of wastewater flow rate observations versus predictions by ARIMA model (a) and MLPNN model (b).

Diagnostic checking

The selected model ARIMA (5,1,1) was also tested and verified by the residual errors of predicted results in addition to the four model evaluation criteria. The time series of residual errors indicates that some observation trends remain to be captured by the model. The distribution of the residual errors suggests the error distribution is Gaussian distribution. The mean value of the residual error is −0.70, which is very close to zero, and the standard deviation is 95.58. The ACF and PACF of residual errors were also examined. ACF and PACF values decay fairly rapidly to zero, either from above or below, and stay within the confidence interval, which implies the selected ARIMA (5,1,1) model can be used to analyze the influent flow rate at the Barrie WwTF.

Results of the MLPNN model

To develop the MLPNN model, the original dataset was also split into two parts: 66% for training and the remaining 34% for testing. The RMSE, MAPE, R2, and NSE for the training set are 115.85, 3.149, 0.791, and 0.747, respectively. For the testing set, the values are 76.63, 2.636, 0.642, and 0.556, respectively. The NSE of the testing set is greater than 0.5, indicating acceptable model performance (Moriasi et al. 2007). Results of observations and predictions obtained by the MLPNN model are plotted in Figure 6(b). It is clearly shown that the MLPNN has good generalization ability and could capture the flow peaks (Figure 6(b)). The scatter plot of MLPNN results (Figure 7(b)) shows that most of the points are close to the diagonal line, which also indicates the MLPNN could be used to predict the quantity of influent wastewater at Barrie WwTF.

Compared to the MLPNN model, the testing set of ARIMA shows lower RMSE and MAPE values and a higher NSE value, which indicates the ARIMA model has a better forecasting performance than the MLPNN model in this case. However, the four evaluation criteria values of both models are close. For WWTP managers, the selection of a suitable model in practice should be based on the availability and quality of data. For example, an ARIMA model is developed using only the flow data, but it requires the flow data to be stationary after certain transformations. Meanwhile, the MLPNN model has no requirement regarding data stationarity, and exogenous variables can be included for training and prediction purposes.

It is worth mentioning that, in this case study, since the range of the training set is larger than that of the testing set, the model's extrapolation ability was not tested. Considering that predictions made beyond the range of the training set tend to be unreliable, if the testing data contain values that exceed the training data range in future applications, the model's extrapolation ability should be evaluated to provide more reliable and reasonable prediction (Pektas & Cigizoglu 2017; Oyebode & Stretch 2019).

CONCLUSIONS

In order to find an effective and efficient way to predict wastewater inflow rate, an ARIMA model and an MLPNN model were developed in this study. The developed models were applied to predict the wastewater inflow rate of the Barrie WwTF. This work attempts to fill the gaps in the wastewater inflow forecasting research in Canada. Both models provided satisfactory forecasting results, which indicates that both of them could be used to predict the influent flow rate at WWTPs similar to the Barrie WwTF. The ARIMA model showed a better forecasting performance than the MLPNN model, in terms of RMSE, MAPE and NSE. However, it should be noted that each model has its own limitations. In the ARIMA model, the time series must be stationary. If the data is still not stationary after certain transformations, ARIMA cannot be used. In the MLPNN model, its performance heavily relies on the quantity and quality of input data. It is not suitable when reliable input data are not easily accessible. Researchers and practitioners should carefully consider these limitations when choosing an appropriate model for wastewater flow forecasting.

ACKNOWLEDGEMENTS

This study was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) and Southern Ontario Water Consortium (SOWC). The authors would also like to acknowledge the City of Barrie for providing the data.

REFERENCES

REFERENCES
Al-Asheh
S.
,
Mjalli
F. S.
&
Alfadala
H. E.
2007
Forecasting influent-effluent wastewater treatment plant using time series analysis and artificial neural network techniques
.
Chemical Product and Process Modeling
2
(
3
).
https://doi.org/10.2202/1934-2659.1063
.
Ashley
R. M.
,
Hvitved-Jacobsen
T.
&
Bertrand-Krajewski
J. L.
1999
Quo vadis sewer process modelling?
Water Science and Technology
39
(
9
),
9
22
.
Benedetti
L.
,
Langeveld
J.
,
Comeau
A.
,
Corominas
L.
,
Daigger
G.
,
Martin
C.
,
Mikkelsen
P. S.
,
Vezzaro
L.
,
Weijers
S.
&
Vanrolleghem
P. A.
2013
Modelling and monitoring of integrated urban wastewater systems: review on status and perspectives
.
Water Science and Technology
68
(
6
),
1203
1215
.
Bennett
N. D.
,
Croke
B. F. W.
,
Guariso
G.
,
Guillaume
J. H. A.
,
Hamilton
S. H.
,
Jakeman
A. J.
,
Marsili-Libelli
S.
,
Newham
L. T. H.
,
Norton
J. P.
,
Perrin
C.
,
Pierce
S. A.
,
Robson
B.
,
Seppelt
R.
,
Voinov
A. A.
,
Fath
B. D.
&
Andreassian
V.
2013
Characterising performance of environmental models
.
Environmental Modelling and Software
40
,
1
20
.
Box
G. E. P.
&
Jenkins
G. M.
1976
Time Series Analysis: Forecasting and Control
.
Holden-Day
,
San Francisco, CA, USA
.
Buchberger
S. G.
&
Wells
G. J.
1996
Intensity, duration, and frequency of residential water demands
.
Journal of Water Resources Planning and Management
122
(
1
),
11
19
.
Campisano
A.
,
Cabot Ple
J.
,
Muschalla
D.
,
Pleau
M.
&
Vanrolleghem
P. A.
2013
Potential and limitations of modern equipment for real time control of urban wastewater systems
.
Urban Water Journal
10
(
5
),
300
311
.
Contreras
J.
,
Espínola
R.
,
Nogales
F. J.
&
Conejo
A. J.
2003
ARIMA models to predict next-day electricity prices
.
IEEE Transactions on Power Systems
18
(
3
),
1014
1020
.
Djebbar
Y.
&
Kadota
P. T.
1998
Estimating sanitary flows using neural networks
.
Water Science and Technology
38
(
10
),
215
222
.
Fernandez
F. J.
,
Seco
A.
,
Ferrer
J.
&
Rodrigo
M. A.
2009
Use of neurofuzzy networks to improve wastewater flow-rate forecasting
.
Environmental Modelling and Software
24
(
6
),
686
693
.
Fletcher
D.
&
Goss
E.
1993
Forecasting with neural networks: an application using bankruptcy data
.
Information and Management
24
(
3
),
159
167
.
Hellerstein
J. M.
2008
Quantitative Data Cleaning for Large Databases
.
United Nations Economic Commission for Europe
. .
Jain
A. K.
,
Mao
J.
&
Mohiuddin
K. M.
1996
Artificial neural network: a tutorial
.
Computer
29
(
3
),
31
44
.
Junninen
H.
,
Niska
H.
,
Tuppurainen
K.
,
Ruuskanen
J.
&
Kolehmainen
M.
2004
Methods for imputation of missing values in air quality data sets
.
Atmospheric Environment
38
(
18
),
2895
2907
.
Kim
J. R.
,
Ko
J. H.
,
Im
J. H.
,
Lee
S. H.
,
Kim
S. H.
,
Kim
C. W.
&
Park
T. J.
2006
Forecasting influent flow rate and composition with occasional data for supervisory management system by time series model
.
Water Science and Technology
53
(
4–5
),
185
192
.
Kim
M.
,
Kim
Y.
,
Kim
H.
,
Piao
W.
&
Kim
C.
2016
Evaluation of the k-nearest neighbor method for forecasting the influent characteristics of wastewater treatment plant
.
Frontiers of Environmental Science and Engineering
10
(
2
),
299
310
.
Lindqvist
J.
,
Wik
T.
,
Lumley
D.
&
Äijälä
G.
2005
Influent load prediction using low order adaptive modelling
. In
2nd IWA Conference on Instrumentation, Control and Automation
,
Busan, South Korea
.
Little
T. D.
, (ed.)
2013
The Oxford Handbook of Quantitative Methods, Volume 2, Statistical Analysis
.
Oxford University Press
,
Oxford, UK
, pp.
458
460
.
Machiwal
D.
&
Jha
M. K.
2012
Hydrologic Time Series Analysis: Theory and Practice
.
Springer
,
Dordrecht
,
The Netherlands
.
McCulloch
W. S.
&
Pitts
W.
1943
A logical calculus of the ideas immanent in nervous activity
.
The Bulletin of Mathematical Biophysics
5
(
4
),
115
133
.
Moriasi
D. N.
,
Arnold
J. G.
,
Van Liew
M. W.
,
Bingner
R. L.
,
Harmel
R. D.
&
Veith
T. L.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Transactions of the ASABE
50
(
3
),
885
900
.
Nielsen
H.
,
Engelbrecht
J.
,
Brunak
S.
&
Heijne
G. V.
1997
A neural network method for identification of prokaryotic and eukaryotic signal peptides and prediction of their cleavage sites
.
International Journal of Neural Systems
8
(
5–6
),
581
599
.
Novotny
V.
,
Jones
H.
,
Feng
X.
&
Capodaglio
A.
1991
Time series analysis models of activated sludge plants
.
Water Science and Technology
23
(
4–6
),
1107
1116
.
Ordieres
J. B.
,
Vergara
E. P.
,
Capuz
R. S.
&
Salazar
R. E.
2005
Neural network prediction model for fine particulate matter (PM2.5) on the US–Mexico border in El Paso (Texas) and Ciudad Juárez (Chihuahua)
.
Environmental Modelling and Software
20
(
5
),
547
559
.
Pektas
A. O.
&
Cigizoglu
H. K.
2017
Investigating the extrapolation performance of neural network models in suspended sediment data
.
Hydrological Sciences Journal
62
(
10
),
1694
1703
.
Rafiq
M. Y.
,
Bugmann
G.
&
Easterbrook
D. J.
2001
Neural network design for engineering applications
.
Computers and Structures
79
(
17
),
1541
1552
.
Schmitt
T. G.
,
Thomas
M.
&
Ettrich
N.
2004
Analysis and modeling of flooding in urban drainage systems
.
Journal of Hydrology
299
(
3–4
),
300
311
.
Tan
P. C.
,
Berger
C. S.
,
Dabke
K. P.
&
Mein
R. G.
1991
Recursive identification and adaptive prediction of wastewater flows
.
Automatica
27
(
5
),
761
768
.
Verma
A.
,
Wei
X.
&
Kusiak
A.
2013
Predicting the total suspended solids in wastewater: a data-mining approach
.
Engineering Applications of Artificial Intelligence
26
(
4
),
1366
1372
.
Wei
X.
,
Kusiak
A.
&
Sadat
H. R.
2013
Prediction of influent flow rate: data-mining approach
.
Journal of Energy Engineering
139
(
2
),
118
123
.