## 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

*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: 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*).

*q*, for which , is the order of the MA process. Thus, the MA process of order

*q*can be denoted by MA(

*q*).

*d*. A process which has

*d*th order differencing transformations is called an integrated process of order

*d*, and it can be denoted as I(

*d*).

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.

#### 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 (R^{2}), 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. R^{2} 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 R^{2} 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 km^{2} 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).

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 m^{3}/day and the design peak flow is 156,000 m^{3}/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

*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.4826

*x*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. where

*k*

_{i}are samples, ,

*n*is the total number of samples, and

*C*= median {

*k*

_{i}}.

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).

C + 4.4478MADs (m^{3}/hr) | C − 4.4478MADs (m^{3}/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% | 2 | 0.02% | 0.59% |

Summer | 3,528.44 | 1,258.62 | 8,832 | 30 | 0.34% | 9 | 0.10% | 0.44% |

Fall | 3,733.94 | 835.06 | 8,736 | 7 | 0.08% | 0 | 0.00% | 0.08% |

Winter | 3,843.89 | 1,073.88 | 8,736 | 0 | 0.00% | 0 | 0.00% | 0.00% |

Year | – | – | 35,132 | 87 | 0.25% | 11 | 0.03% | 0.28% |

C + 4.4478MADs (m^{3}/hr) | C − 4.4478MADs (m^{3}/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% | 2 | 0.02% | 0.59% |

Summer | 3,528.44 | 1,258.62 | 8,832 | 30 | 0.34% | 9 | 0.10% | 0.44% |

Fall | 3,733.94 | 835.06 | 8,736 | 7 | 0.08% | 0 | 0.00% | 0.08% |

Winter | 3,843.89 | 1,073.88 | 8,736 | 0 | 0.00% | 0 | 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 m^{3}/hr to 4,206 m^{3}/hr, with an average rate of 2,335 m^{3}/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 m^{3}/hr in the spring due to snow melting, whereas the lowest inflow recorded was 1,056 m^{3}/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.

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 m^{3}/hr) in March is much broader than that in April and May (1,297–3,940 and 1,180–3,565 m^{3}/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 m^{3}/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 m^{3}/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 m^{3}/hr, except for Friday with a maximum of 2,739 m^{3}/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 m^{3}/hr, then began to increase to reach a maximum of 3,116 m^{3}/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.

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.

ARIMA (p, d, q) | Training | Testing | ||||||
---|---|---|---|---|---|---|---|---|

RMSE | MAPE (%) | R^{2} | NSE | RMSE | MAPE (%) | R^{2} | 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 (%) | R^{2} | NSE | RMSE | MAPE (%) | R^{2} | 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 R^{2}, and highest NSE values are marked bold.

#### Model validation

As shown in Table 2, the RMSE, MAPE (%), R^{2}, 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.

#### 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, R^{2}, 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

*United Nations Economic Commission for Europe*