## Abstract

The management of wastewater treatment plant (WWTP) and the assessment of uncertainty in its design are crucial from an environmental engineering perspective. One of the key mechanisms in WWTP operation is activated sludge, which is related to the biological oxygen demand (BOD) parameter. In the modeling of BOD, the conventional approach utilizing ordinary differential equations (ODEs) fails to incorporate the stochastic nature of this parameter, leading to a considerable level of uncertainty in the design of WWTP. To address this issue, this study proposes a stochastic model that utilizes stochastic differential equations (SDEs) instead of ODE to simulate BOD activities of microorganisms and wastewater flow rate (*Q*). The SDEs and integral It̂o are solved using the Euler–Maruyama method for a period of 15 sequential days and the timespan of 2019–2020 for a WWTP in Tabriz City. SDE improves the design of WWTP facilities by identifying uncertainties and enhancing reliability. This, in turn, increases the reliability of the technical structures within the WWTP and improves the performance of its biological system. By considering the randomness of the problem, the proposed method significantly improves the results, with an enhancement of 11.47 and 10.11% for the BOD and *Q* models, respectively.

## HIGHLIGHTS

This study suggests using stochastic differential equations (SDEs) and the Euler–Maruyama method instead of ordinary differential equations (ODEs) to model BOD and

*Q*in WWTP design, addressing uncertainty.A 15-day simulation was conducted to evaluate the variations of BOD and

*Q*pollution parameters in the biological section of WWTP for the first time.

## ABBREVIATIONS

- WWTP
Wastewater treatment plant

- BOD (or BOD
_{5}) Total 5-day biochemical oxygen demand

- COD
Total chemical oxygen demand

- MLSS
Mixed liquid suspended solids

- SDE
Stochastic differential equations

- Q
Flow rate of wastewater

- ODE
Ordinary differential equations

- inf
Sampled at the input of the biological reactor

- eff
Sampled at the effluent of the biological reactor

- PH
Potential of hydrogen

- SS
Suspended solids

- TDS
Total dissolved solids

- TSS
Total suspended solids

- VSS
Volatile suspended solids

- EC
Electrical conductivity

## SYMBOLS

Reactor volume

*T*Air temperature

Biomass concentration

The final concentration of biomass

Return activated sludge flow rate

Activated sludge effluent from the WWTP system

The average flow rate

Endogenous decay coefficient

Effluent substrate concentration

Influent substrate concentration

The average substrate concentration

The cell efficiency factor

Mean cell residence time (MCRT)

*r*Sampled at returning sludge

Wiener process

*f*The drift coefficient

*g*The diffusion coefficient

*a*,*b*and*c*Numerical constant coefficients

*S*_{t}An explicit answer of the Black–Scholes equation

*μ*Mean of data

*σ*^{2}Variance of data

## INTRODUCTION

Managing water resources for human societies requires the utilization of urban wastewater as a part of the sustainable water cycle. Wastewater reuse has been acknowledged as a fundamental strategy in safeguarding water resources and providing solutions for global water challenges (Wester *et al.* 2015; Gul *et al.* 2021). To achieve this objective, it is crucial to design close-to-reality wastewater treatment plants (WWTP) acknowledging that the operation and performance of any WWTP are significantly defined by the quality of both the untreated and treated effluents (Ahmadi *et al.* 2017; Nourani *et al.* 2018). Therefore, it is necessary to establish a direct connection (by establishing a more precise scientific link) with the treatment plant's output for industrial, agricultural, public, or sanitary purposes in order to meet environmental standards and indicators (Cresson & Sonner 2018). This will help assess uncertainty regarding physical and chemical parameters and enhance productivity. However, modeling the performance of a WWTP is challenging due to the complexity of treatment processes (Hanbay *et al.* 2008; Nourani *et al.* 2023a, 2023b). In order to enhance the design quality of treated wastewater, it is crucial to monitor and stabilize the physical variables associated with the proper functioning of activated sludge in WWTP. Consequently, applying stochastic mathematical equations to model sewage output, which can be considered a form of pollution control problem, has the potential to minimize uncertainty and enhance the accuracy of the models. Advancing our understanding of wastewater output can contribute significantly to resolving pollution control challenges.

The quality of the effluent of the WWTP undergoes a transformation through the activated sludge process, and this alteration can be assessed by measuring the biological oxygen demand (BOD). Conceptual models with mathematical relationships are commonly used to model BOD in the biological part of WWTP. However, there are two major problems that make modeling difficult when it comes to the mathematical modeling of environmental systems based on the physics of the problem. Firstly, the multitude of variables increases the number of options to be examined and the number of variables to be considered. Secondly, assessing uncertainty becomes challenging due to the intricate connections among the components engaged in stochastic process modeling and the significant level of randomness involved. In order to model BOD and the design of WWTP, ordinary differential equations (ODEs), the equations involving ordinary derivatives of a function, are typically used as mathematical solutions due to their ease of use (e.g., Diehl & Farås 2013; Mauritsson 2013). However, most natural processes are stochastic in nature, and the uncertainty caused by randomness needs to be taken into account when modeling such processes (Nourani *et al.* 2023a, 2023b). The ODE model alone cannot scientifically interpret the true physics of the phenomenon so ignore the uncertainties. Therefore, using the ODE model, uncertainties in the simulated real stochastic systems are disregarded (Compagnoni *et al.* 2023).

Numerical simulations, such as those performed using ODEs, often require accurate measurements or the adjustment of mathematical models based on prototype measurements from material models (Wylie & Streeter 1978). There are three classifications for hydrological/mathematical models: (1) black-box models, (2) conceptual models (sometimes called gray-box models), and (3) physically based models. In the physical approach (sometimes called white-box models or theoretical models), to estimate the output of a real system, it is necessary to pay attention to the inputs and the structure of the system and also its governing physical laws. In such a way, the initial and boundary conditions serve as essential constraints for solving differential equations that govern the modeling. Mathematical models based on a physical approach incorporate uncertainty stemming from simplification and assumptions made during the modeling process. Additionally, they may be susceptible to uncertainty arising from the collected physical data. The simulation of a real system can be achieved by utilizing the numerical model of the governing equations. For instance, in the analysis of activated sludge, the law of conservation of mass is applied to the input and output components of the biological part of a WWTP where the governing equations in WWTP design were of the ODE type (Metcalf 2003; Sin *et al.* 2009; Lin 2014). However, wastewater treatment using an activated sludge process involves a combination of mechanical, chemical, and biological actions, and can be monitored using nonlinear and complete equations with discontinuous solutions (Mauritsson 2013).

To assess uncertainties, investigations in the field of designing and simulating real random phenomena have continued in advanced mathematics on stochastic spaces by handling stochastic differential equations (SDEs) (Dexter *et al.* 2016). SDEs are a better alternative to ODEs for simulating random processes by taking into account uncertainties as well (Kabouris & Georgakakos 1991). An SDE is a differential equation that is presented along with a random variable, such as white noise. In this way, SDEs are used to describe Brownian Motion and Weiner Process. In SDE, one or more of the terms is a stochastic process, resulting in a solution that is also a stochastic process. A stochastic model through SDEs includes random variables that change over time. These variables can be described by a probability distribution. Stochastic modeling helps estimate variations in a system's state and output (Ouldali *et al.* 1989). SDEs are widely used in complex modeling, such as in the physical modeling of WWTP (e.g., Beck 1987; Browning *et al.* 2020). For instance, Harris (1977) investigated SDEs using Euler equations for WWTP with the activated sludge system, considering daily data over 2 sequential days. His research was innovative, and a stabilizing pool at the influent of the treatment plant was considered to apply a constant flow rate and eliminate stochastic values. To simplify the simulation, most of the parameters and numerical values were considered constant. However, the results showed a significant difference between the calculated and observed BOD. Given the alterations of 100 m^{3}/day for the inflow rate and 240 mg/L for BOD, in the context of a 2-day scale investigation and the presence of the scale effect in these computations, the error in the calculations was also noteworthy.

The current study utilized SDEs to analyze Brownian motion in simulating WWTP in Tabriz City. To solve the stochastic processes of BOD and *Q* by SDE solution, the integral and Euler–Maruyama method were used. The integral analyses a stochastic process with respect to Brownian motion within a bounded time interval [0, *t*] (Erfanian *et al.* 2016). Optimal coefficients of SDE were extracted by comparing computational and observational data over 2 weeks (the longest modeling run for BOD and *Q* variables of WWTP lasted 2 consecutive days, based on prior research by Harris (1977)). The application of this method in the assessment of uncertainty in the activated sludge component is considered an innovation to the present research. A Wiener model using SDE was presented to check real data and measure pollutant biological parameters for short-term changes in BOD and *Q*, and to assess uncertainties in the biological reactor section.

## MATERIALS AND METHODS

### Study site and data

The study examined the assessment of uncertainty in simulating BOD and *Q* agents in WWTP design, with a focus on SDE solutions. Simulations ran from 2019 to 2020 for 15 days (for the first time), for the BOD pollution index and the amount of *Q* changes in WWTP's biological part in Tabriz City with an activated sludge system.

^{3}s

^{−1}, but the actual average discharge is 1.15 m

^{3}s

^{−1}. Records indicate that the maximum discharge on non-rainy days reaches 2.5 m

^{3}s

^{−1}, while on rainy days, it can go up to 3.8 m

^{3}s

^{−1}. The influent concentrations of BOD and total chemical oxygen demand (COD) variables are 250 and 320 mg l

^{−1}, respectively. The Tabriz City WWTP consists of physical, biological, and chemical treatment stages. The incoming sewage carries suspended materials weighing 32,400 and 45,360 kg day

^{−1}, respectively. Table 1 presents various statistical characteristics of WWTP's pollution indexes. Also, as target parameters in SDE modeling, Figure 2 illustrates the daily fluctuations of

*Q*

_{inf}and

*Q*

_{eff}to WWTP over a 1-year period (2019–2020).

Variable . | Measurement stage . | Unit . | Minimum . | Maximum . | Mean . | Std . |
---|---|---|---|---|---|---|

Influent to WWTP | – | 7.29 | 8.61 | 8.05 | 0.19 | |

Influent to WWTP | mg L^{−1} | 220.00 | 315.00 | 290.12 | 13.05 | |

Influent to WWTP | mg L^{−1} | 362.17 | 516.67 | 479.16 | 22.01 | |

Influent to WWTP | mg L^{−1} | 0.00 | 1.50 | 0.61 | 0.151 | |

Influent to WWTP | mg L^{−1} | 234.00 | 336.33 | 296.40 | 14.64 | |

Return activated sludge | mg L^{−1} | 3,048.67 | 3,892.00 | 3,271.21 | 89.83 | |

Influent to WWTP | mg L^{−1} | 140.50 | 235.33 | 191.13 | 15.12 | |

Influent to WWTP | 1,067.17 | 1,366.83 | 1,202.01 | 51.19 | ||

Influent to WWTP | °C | 14.30 | 26.78 | 21.04 | 2.98 | |

Air temperature | °C | −9.08 | 31.42 | 12.42 | 10.84 | |

Effluent from WWTP | mg L^{−1} | 17.33 | 26.33 | 22.72 | 1.28 | |

Effluent from primary settling | mg L^{−1} | 140.00 | 197.50 | 181.37 | 8.48 | |

Secondary sedimentation effluent | mg L^{−1} | 17.50 | 27.33 | 23.66 | 1.35 | |

Effluent from primary settling | mg L^{−1} | 235.33 | 320.50 | 297.25 | 13.33 | |

Effluent from WWTP | mg L^{−1} | 30.00 | 40.83 | 37.98 | 1.87 | |

Effluent from WWTP | 779.83 | 1,188.50 | 936.34 | 64.76 | ||

Effluent from WWTP | mg L^{−1} | 6.33 | 13.13 | 9.68 | 1.16 | |

Effluent from WWTP | °C | 13.97 | 27.73 | 21.05 | 3.34 | |

Effluent from primary settling | mg L^{−1} | 96.67 | 132.33 | 119.46 | 5.55 | |

Secondary sedimentation effluent | mg L^{−1} | 14.80 | 26.45 | 22.30 | 2.09 | |

Effluent from WWTP | mg L^{−1} | 14.47 | 25.98 | 21.92 | 2.043 | |

Effluent from primary settling | mg L^{−1} | 235.33 | 320.50 | 297.25 | 13.33 | |

Effluent from primary settling | – | 7.09 | 8.16 | 7.70 | 0.128 | |

Return activated sludge | – | 6.72 | 7.95 | 7.43 | 0.156 | |

Active sludge | mg L^{−1} | 1,138.00 | 2,870.00 | 2,273.48 | 147.50 | |

Active sludge | μSiemens cm^{−1} | 824.00 | 2,320.00 | 1,550.23 | 280.76 |

Variable . | Measurement stage . | Unit . | Minimum . | Maximum . | Mean . | Std . |
---|---|---|---|---|---|---|

Influent to WWTP | – | 7.29 | 8.61 | 8.05 | 0.19 | |

Influent to WWTP | mg L^{−1} | 220.00 | 315.00 | 290.12 | 13.05 | |

Influent to WWTP | mg L^{−1} | 362.17 | 516.67 | 479.16 | 22.01 | |

Influent to WWTP | mg L^{−1} | 0.00 | 1.50 | 0.61 | 0.151 | |

Influent to WWTP | mg L^{−1} | 234.00 | 336.33 | 296.40 | 14.64 | |

Return activated sludge | mg L^{−1} | 3,048.67 | 3,892.00 | 3,271.21 | 89.83 | |

Influent to WWTP | mg L^{−1} | 140.50 | 235.33 | 191.13 | 15.12 | |

Influent to WWTP | 1,067.17 | 1,366.83 | 1,202.01 | 51.19 | ||

Influent to WWTP | °C | 14.30 | 26.78 | 21.04 | 2.98 | |

Air temperature | °C | −9.08 | 31.42 | 12.42 | 10.84 | |

Effluent from WWTP | mg L^{−1} | 17.33 | 26.33 | 22.72 | 1.28 | |

Effluent from primary settling | mg L^{−1} | 140.00 | 197.50 | 181.37 | 8.48 | |

Secondary sedimentation effluent | mg L^{−1} | 17.50 | 27.33 | 23.66 | 1.35 | |

Effluent from primary settling | mg L^{−1} | 235.33 | 320.50 | 297.25 | 13.33 | |

Effluent from WWTP | mg L^{−1} | 30.00 | 40.83 | 37.98 | 1.87 | |

Effluent from WWTP | 779.83 | 1,188.50 | 936.34 | 64.76 | ||

Effluent from WWTP | mg L^{−1} | 6.33 | 13.13 | 9.68 | 1.16 | |

Effluent from WWTP | °C | 13.97 | 27.73 | 21.05 | 3.34 | |

Effluent from primary settling | mg L^{−1} | 96.67 | 132.33 | 119.46 | 5.55 | |

Secondary sedimentation effluent | mg L^{−1} | 14.80 | 26.45 | 22.30 | 2.09 | |

Effluent from WWTP | mg L^{−1} | 14.47 | 25.98 | 21.92 | 2.043 | |

Effluent from primary settling | mg L^{−1} | 235.33 | 320.50 | 297.25 | 13.33 | |

Effluent from primary settling | – | 7.09 | 8.16 | 7.70 | 0.128 | |

Return activated sludge | – | 6.72 | 7.95 | 7.43 | 0.156 | |

Active sludge | mg L^{−1} | 1,138.00 | 2,870.00 | 2,273.48 | 147.50 | |

Active sludge | μSiemens cm^{−1} | 824.00 | 2,320.00 | 1,550.23 | 280.76 |

### Stochastic differential equation

The probability space associated with a stochastic phenomenon is represented by a continuous path Wiener process within a specific time interval , characterized by independent variables. As a result of the presence of Wiener processes in the corresponding stochastic integral, the solution to these equations is not derivative-based concerning time (*t*), and typically, such equations are expressed in terms of stochastic integrals . The stochastic integral is a concept within stochastic calculus that expands upon the Riemann–Stieltjes integral (which extends the Riemann integral to integrate functions with respect to another function) to encompass stochastic processes. This integral is of great significance in the study of SDEs and finds utility in various disciplines, such as finance, economics, and physics. By accounting for the stochastic nature and volatility of the underlying process, the stochastic integral is defined with respect to a continuous stochastic process. It is employed to solve SDEs and analyze the behavior of stochastic processes (Fima 1999).

*et al.*(2016):where in a

*d*-dimensional field, let be a set of variables,

*f*be the drift coefficient,

*g*be the diffusion coefficient of the (

*d*×

*m*)-dimensional matrix, and be the

*m*-dimensional Wiener process. To define the integral with a certain time step and the condition a part of the range , and are considered. In this case, the integral can be defined as shown in the following equation:

*g*is a measurable function and the integral of is finite, then the limit of Equation (2) exists as a result of the integral. The Euler–Maruyama method is used for the numerical solution of SDE in this study. Geometric Brownian motion is often used in technical analysis, which includes positive and negative variables (Erfanian

*et al.*2016). When modeling the parameters of the treatment plant, the activated sludge section, and the aeration lagoons, it is important to consider that the system is progressive and the measured variables will not be negative. Therefore, Brownian motion with log-normal distribution, which yields only non-negative values, was used. The changes in in this system follow the recursive process of the following equation.where

*a*,

*b*, and

*c*are numerical constant coefficients that follow the Hull and White risk natural stochastic volatility model (Higham 2001). represents changes in Brownian motion, and is an explicit answer of the Black–Scholes equation, which is a more general form of Equation (1). The value of can be determined using the following equation:

Symbols *μ* and *σ*^{2} represent the mean and variance of data, respectively. *σ* denotes standard deviation and *W _{t}* represents Wiener probability process at each time step of

*δt*.

#### SDE for activated sludge process: the concept and algorithm

*et al*. 2022). Therefore, in this study, SDE was utilized to mitigate the impact of uncertainties related to the activated sludge stage and the biological parameters of the effluent or influent system. The activated sludge process is formed in the secondary sedimentation section (Metcalf 2003). The variables involved in the formation of the activated sludge process are depicted in Figure 4. Within this cycle,

*Q*represents the secondary influent flow rate (m

^{3}day

^{−1}),

*Q*

_{R}denotes the return activated sludge flow rate (m

^{3}day

^{−1}) with concentration , and

*Q*

_{w}signifies the activated sludge effluent from the WWTP system (m

^{3}day

^{−1}). ∀ represents the volume of the tank or reactor (m

^{3}),

*X*represents the biomass concentration (mg L

^{−1}), represents the final concentration of biomass (mg L

^{−1}),

*S*

_{0}represents the influent substrate concentration at time

*t*

_{0}, which is equivalent to BOD

_{inf}, and

*S*

_{e}represents the effluent substrate concentration at time

*t*, which is equivalent to (mg L

^{−1}) as provided by the environmental organization. The primary objective of this unit is to separate a large volume of mixed liquor suspended solids (MLSS) from the aeration tank and produce a clear and stable stream with low total suspended solids (TSS) and BOD concentration. This cycle is illustrated schematically in Figure 4.

#### SDE for activated sludge process: math relationships development

*y*represents the cell efficiency factor, which depends on substrate characteristics, activated sludge type, and temperature. represents the endogenous decay coefficient (

*T*

^{−1}). In the second step, the treatment unit's capacity is determined using the following equation:where

*t*is the hydraulic retention time (HRT), and ∀ is the tank volume (m

^{3}). represents the mean cell residence time (MCRT), a crucial design factor for activities of microorganisms in the biological reactor, obtainable from the following equation.the value of must exceed the duration of the highest bacterial population growth in the system in order to avoid the release of biomass or biological clots with the effluent. Assuming a constant flow rate, the equation of continuity in the biological reactor section yields the following equation:the term (

*Q*+

*Q*

_{R})

*X*represents the influential biological flocculation, while the term (

*Q*

_{w}+

*Q*

_{R})

*X*

_{u}represents the effluent biological flocculation (mg L

^{−1}). The inputs and outputs for these equations are listed in Table 2.

Variable . | Inputs . | Outputs . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

. | . | . | χ_{u}
. | χ
. | . | . | . | . | . | |

Units | m^{3} | day^{−1} | – | mg L^{−1} | mg L^{−1} | mg L^{−1} | m^{3} day^{−1} | m^{3} day^{−1} | m^{3} day^{−1} | mg L^{−1} |

Variable . | Inputs . | Outputs . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

. | . | . | χ_{u}
. | χ
. | . | . | . | . | . | |

Units | m^{3} | day^{−1} | – | mg L^{−1} | mg L^{−1} | mg L^{−1} | m^{3} day^{−1} | m^{3} day^{−1} | m^{3} day^{−1} | mg L^{−1} |

Hence, the initial objective of the study was to derive the Wiener form of Equation (5) for *Q*_{eff} and BOD_{eff} in the first phase, and subsequently transform them into stochastic form. It means that the equations pertaining to activated sludge were derived, taking into account the uncertainties associated with the stochastic mode and the Euler–Maruyama form. In the second phase, the relevant analytical codes were randomly executed over a period of 15 days, spanning from 2019 to 2020, and the outcomes were documented. The data used for the model were selected from a 1-year period. Since it was possible to analyze the data using SDE, it was divided into periods of up to 15 days. It was not possible to increase the time period for more than 15 days, in order to apply the divergence in the calculations and account for the scale effect. The calculations are divided into 2-week periods for 1 year. After various modeling and optimal value extraction, the best model was chosen. The calculations incorporated variations in instantaneous time steps and daily data collection, owing to the inherent uncertainty in the physics of the phenomenon. The stochastic variations can be analyzed from three different perspectives:

*Q*_{inf}and BOD_{inf}uncertainties, and their influence on*X,*impacted by climate change, precipitation, and flow volume, can decrease wastewater BOD. Consumption patterns during seasons can also affect these parameters.The uncertainty surrounding the parameters that govern the behavior of microorganisms, specifically

*k*_{d}and*y*, as well as the impact of past data on these values, is a matter of concern. The daily fluctuations in the birth and death rates of the population have a direct bearing on the aforementioned parameters, and their investigation is imperative in light of the uncertainties involved.Flow rates fluctuate in certain hours of the day to their maximum and minimum values over a period of 1 week to 1 month, but are treated as constant in calculations and significantly affect activated sludge parameters.

*dS*. The following equation connects the ratio of to .

In this context, and represent mass entering the system, while represents mass exiting the system. is zero at the influent stage of the WWTP. The effluent mass from the system is denoted as , which has a negative connotation. The parameter is associated with the mass of the reactor. In Equations (12) and (13), the values of and are unknown. represents the quantity of sludge, while represents the quantity of . is a limited value determined according to environmental organization standards.

*X*

_{u}is influenced by sludge concentration. The final step is stochastic equations based on

*Q*and

*S*

_{e}(or BOD). A sequence of variables

*Q*and produces data to derive Equations (14) and (15) representing stochastic equations for flow rate and sludge concentration. The Euler–Maruyama technique with geometric Brownian motion and a log-normal distribution is used to derive these equations.

is the average flow rate, is the average BOD based on observational data, and are specific values at a given moment, are stochastic independent variables in the form of following a Brownian motion with a log-normal distribution, are numerical constants determined through trial and error, and *σ* is the standard deviation of the observational data.

*Q*variables using the proposed model (Scott

*et al.*2003).where

*x*represents the ODE design of the WWTP (Metcalf 2003),

*y*stands for the average SDE simulation result, and

*z*denotes the average observational data. The disparity between

*r*1 and

*r*2, i.e., , determines the improvement percentage.

## RESULTS

### SDE modeling results for flow rate (*Q*)

The assessment of uncertainty took into account the fluctuations in influent flow caused by varying consumption patterns. To prevent turbulent flows during floods or rains, a flow stabilization criterion was implemented, ensuring a specific range of fluctuations. This characteristic results in the imposition of a specific level of influent sewage load on the treatment facility. However, uncertainties arising from surface runoff or other factors could not be measured at specific time intervals for excessive design loads. The coefficients in Equation (14) were determined to match the observational data over a 2-week computational time series. The Wiener process (*W _{i}*) is employed to model uncertainties, with the vector function

*a*and matrix function

_{i}*c*being influenced by several factors such as the flow rate and volume of sludge, as well as the overall kinetic and stoichiometric biological constants found in wastewater. The study focuses on a real model of a WWTP and Table 3 presents the modeling results, comparing the optimal coefficients of stochastic equations with the observational data over a period of 15 days.

_{i}Day . | Observation data (m^{3} day^{−}^{1})
. | Random number . | . | . | dw
. _{i} | Prediction data (m^{3} day^{−1})
. | Constant coefficients . |
---|---|---|---|---|---|---|---|

1 | 142,719 | 0.188 | 57.87 | 377.78 | 0.26 | 142,719.00 | a_{1} = 0.0035 |

2 | 142,691 | 0.910 | −85.57 | 377.78 | −0.42 | 1,427,760.87 | c_{1} = 0.0028 |

3 | 142,653 | 0.484 | 89.70 | 377.78 | 0.41 | 142,689.29 | dt= 1 × 10^{−3} [s] |

4 | 142,663 | 0.532 | 10.36 | 377.80 | 0.04 | 142,779.00 | |

5 | 142,677 | 0.711 | −66.98 | 377.84 | −0.32 | 142,789.35 | |

6 | 142,688 | 0.694 | −15.87 | 377.83 | −0.08 | 142,722.38 | |

7 | 142,639 | 0.119 | −143.50 | 377.87 | −0.68 | 142,706.51 | |

8 | 142,649 | 0.037 | 89.52 | 377.78 | 0.41 | 142,563.00 | |

9 | 142,604 | 0.681 | −26.27 | 377.80 | −0.13 | 142,652.53 | |

10 | 142,673 | 0.322 | 63.94 | 377.76 | 0.29 | 142,626.25 | |

11 | 142,719 | 0.367 | 27.88 | 377.65 | 0.12 | 142,690.19 | |

12 | 142,667 | 0.931 | −22.33 | 377.62 | −0.11 | 142,718.08 | |

13 | 142,555 | 0.390 | −2.83 | 377.64 | −0.02 | 142,695.74 | |

14 | 142,694 | 0.908 | −2.83 | 377.68 | 0.28 | 142,692.91 | |

15 | 142,611 | 0.019 | −10.33 | 377.75 | −0.06 | 142,754.81 | |

16 | 142,673 | 0.261 | 65.27 | 377.79 | 0.29 | 142,744.48 | |

α = 42.98 _{Q} |

Day . | Observation data (m^{3} day^{−}^{1})
. | Random number . | . | . | dw
. _{i} | Prediction data (m^{3} day^{−1})
. | Constant coefficients . |
---|---|---|---|---|---|---|---|

1 | 142,719 | 0.188 | 57.87 | 377.78 | 0.26 | 142,719.00 | a_{1} = 0.0035 |

2 | 142,691 | 0.910 | −85.57 | 377.78 | −0.42 | 1,427,760.87 | c_{1} = 0.0028 |

3 | 142,653 | 0.484 | 89.70 | 377.78 | 0.41 | 142,689.29 | dt= 1 × 10^{−3} [s] |

4 | 142,663 | 0.532 | 10.36 | 377.80 | 0.04 | 142,779.00 | |

5 | 142,677 | 0.711 | −66.98 | 377.84 | −0.32 | 142,789.35 | |

6 | 142,688 | 0.694 | −15.87 | 377.83 | −0.08 | 142,722.38 | |

7 | 142,639 | 0.119 | −143.50 | 377.87 | −0.68 | 142,706.51 | |

8 | 142,649 | 0.037 | 89.52 | 377.78 | 0.41 | 142,563.00 | |

9 | 142,604 | 0.681 | −26.27 | 377.80 | −0.13 | 142,652.53 | |

10 | 142,673 | 0.322 | 63.94 | 377.76 | 0.29 | 142,626.25 | |

11 | 142,719 | 0.367 | 27.88 | 377.65 | 0.12 | 142,690.19 | |

12 | 142,667 | 0.931 | −22.33 | 377.62 | −0.11 | 142,718.08 | |

13 | 142,555 | 0.390 | −2.83 | 377.64 | −0.02 | 142,695.74 | |

14 | 142,694 | 0.908 | −2.83 | 377.68 | 0.28 | 142,692.91 | |

15 | 142,611 | 0.019 | −10.33 | 377.75 | −0.06 | 142,754.81 | |

16 | 142,673 | 0.261 | 65.27 | 377.79 | 0.29 | 142,744.48 | |

α = 42.98 _{Q} |

*Q*or BOD have a wide range of changes, and the scale effect has a significant impact on this analysis as well. As a result, abnormal jumps and extreme fluctuations were observed in the results within the 15-day time period, with considerable intervals between them. It means that the convergence of answers obtained for long time intervals led to the calculations were time-consuming. In order to determine the optimal coefficients of the SDEs, a significant amount of trial and error has been employed. These limitations are important to note when considering the method presented in this research.

*Q*

_{eff}for each day, encompassing both the maximum and minimum values, are illustrated in Figure 7 based on the results.

### SDE modeling results for BOD

*t*to

_{i}*t*

_{i}_{+1}. The findings were subsequently presented in Table 4.

Day . | Observation data (mg L^{−1})
. | Random number . | . | . | dw_{i}
. | Prediction data (mg L^{−1})
. | Constant coefficients . |
---|---|---|---|---|---|---|---|

1 | 24 | 0.511 | −0.9375 | 4.8990 | −0.338 | 24.000 | a_{1} = 0.005 |

2 | 23 | 0.307 | −0.5002 | 4.8541 | 0.000 | 20.893 | c_{1} = 0.7 |

3 | 24 | 0.618 | −0.0346 | 4.8059 | 0.110 | 21.014 | dt = 1 × 10^{−3} [s] |

4 | 25 | 0.584 | −0.1654 | 4.8195 | −0.343 | 22.058 | |

5 | 22 | 0.350 | −0.2729 | 4.8307 | 0.267 | 19.161 | |

6 | 23 | 0.577 | −0.8680 | 4.8919 | −0.300 | 21.316 | |

7 | 21 | 0.761 | 0.3462 | 4.7662 | −0.029 | 18.882 | |

8 | 23 | 0.966 | 0.1595 | 4.7857 | 0.677 | 18.776 | |

9 | 22 | 0.481 | 0.7090 | 7.7279 | 0.090 | 23.968 | |

10 | 23 | 0.655 | 0.7498 | 7.7236 | 0.213 | 24.969 | |

11 | 24 | 0.534 | 0.9054 | 4.7071 | −0.020 | 27.238 | |

12 | 25 | 0.315 | 1.0441 | 4.6924 | 0.077 | 27.179 | |

13 | 22 | 0.653 | 0.5534 | 4.7444 | −0.504 | 28.171 | |

14 | 23 | 0.998 | 1.7057 | 4.6213 | −0.199 | 22.653 | |

15 | 21 | 0.394 | 2.1998 | 4.5676 | 0.356 | 20.980 | |

16 | 24 | 0.740 | 1.8310 | 4.6078 | −0.383 | 24.089 | |

Day . | Observation data (mg L^{−1})
. | Random number . | . | . | dw_{i}
. | Prediction data (mg L^{−1})
. | Constant coefficients . |
---|---|---|---|---|---|---|---|

1 | 24 | 0.511 | −0.9375 | 4.8990 | −0.338 | 24.000 | a_{1} = 0.005 |

2 | 23 | 0.307 | −0.5002 | 4.8541 | 0.000 | 20.893 | c_{1} = 0.7 |

3 | 24 | 0.618 | −0.0346 | 4.8059 | 0.110 | 21.014 | dt = 1 × 10^{−3} [s] |

4 | 25 | 0.584 | −0.1654 | 4.8195 | −0.343 | 22.058 | |

5 | 22 | 0.350 | −0.2729 | 4.8307 | 0.267 | 19.161 | |

6 | 23 | 0.577 | −0.8680 | 4.8919 | −0.300 | 21.316 | |

7 | 21 | 0.761 | 0.3462 | 4.7662 | −0.029 | 18.882 | |

8 | 23 | 0.966 | 0.1595 | 4.7857 | 0.677 | 18.776 | |

9 | 22 | 0.481 | 0.7090 | 7.7279 | 0.090 | 23.968 | |

10 | 23 | 0.655 | 0.7498 | 7.7236 | 0.213 | 24.969 | |

11 | 24 | 0.534 | 0.9054 | 4.7071 | −0.020 | 27.238 | |

12 | 25 | 0.315 | 1.0441 | 4.6924 | 0.077 | 27.179 | |

13 | 22 | 0.653 | 0.5534 | 4.7444 | −0.504 | 28.171 | |

14 | 23 | 0.998 | 1.7057 | 4.6213 | −0.199 | 22.653 | |

15 | 21 | 0.394 | 2.1998 | 4.5676 | 0.356 | 20.980 | |

16 | 24 | 0.740 | 1.8310 | 4.6078 | −0.383 | 24.089 | |

_{eff}data for a consecutive period of 15 days, in conjunction with

*Q*

_{eff}, spanning from 2019 to 2020. The outcomes of the best model are depicted in Figure 9, which illustrates a comparison between the computational and observational data. The group of responses for BOD

_{eff}was determined similar to

*Q*

_{eff}calculations. The calculation process involved a time period of 15 days and 10 calculation steps, which were determined through trial and error. The first day's data were a fixed number from observation data. Figure 10 depicts a variation of the BOD

_{eff}during SDE modeling, including the highest, mean, and lowest values obtained for each day.

## DISCUSSION

To design WWTP, using linear analyses overlooks important factors of fluctuations in wastewater phenomena such as critical rainfall events. Insufficient understanding of the uncertainties involved in the processes often leads to significant errors in the calculation factors and implementation volumes of WWTP projects, resulting in increased financial burdens. For instance, failure to accurately calculate the minimum coefficient *Q*_{inf} to design with ODE in critical situations can lead to prolonged hydraulic retention times in sedimentation ponds, aeration, and the time required for settled sludge to reach the nitrification stage. This, in turn, causes an increase in the sedimentation level of compressed particles, resulting in a decline in the quality of treated wastewater. In this study, the foundation of the work was established upon models for stochastic processes, which have the potential to make a significant contribution to the development of models based on linear first-order ODEs, utilized for the analysis of equations in the activated sludge section of the WWTP. To achieve this objective, the stochastic form of the equations for two output parameters, namely *Q*_{eff} and BOD_{eff}, was extracted from the series of linear equations of the activated sludge section. Afterwards, the constant coefficients (*a _{i}* and

*c*) were modified for both models. The optimal values for

_{i}*a*and

_{i}*c*for

_{i}*Q*

_{eff}were 0.0035 and 0.0028, respectively, while for BOD

_{eff}, they were obtained as 0.005 and 0.7, respectively. To ensure maximum convergence in results,

*δt*was set to 0.001. However, one of the limitations of this method is its inability to accurately compute for long-term time scales. The presence of computation errors in each time step leads to the accumulation of unacceptable errors in modeling. Nevertheless, this method demonstrated a significant ability to monitor the fluctuations and activities of microorganisms and control the changes in short time scales. It is noteworthy that the used method can handle uncertainty in the activated sludge process for short time scales, while the results for long scales may not be reliable.

Due to the absence of analytical solutions, numerical methods have been necessary for the computational simulation of most SDEs, especially nonlinear ones. In the modeling process, the results at each step, *t _{i}* compared to

*t*, exhibit variations due to the presence of Brownian motion and

_{i+1}*W*values. The objective is to achieve a numerical approximation that minimizes the error between observational and computational data, thereby obtaining optimal constant coefficients through such a solution. In this regard, analytical results demonstrate the convergence of the Euler–Maruyama method for small-time step sizes. The typical fluctuations observed in evaluating SDEs are practically absent in the activated sludge process, which is considered a macro physical phenomenon that eliminates Wiener process noises. The models developed based on SDE theory, as integrated noisy models, serve two fundamental purposes. Firstly, it is crucial to ascertain that the solution process exhibits a Markovian characteristic, wherein the value of

_{i}*S*

_{i+}_{1}is contingent upon

*S*. Secondly, the set of obtained results must belong to a broad class of martingales, guaranteeing that the solution process exhibits the characteristics of Brownian motion or Wiener process.

_{i}A comprehensive comparison between the final performance of ODE and SDE analysis with observational data is presented in Table 5. In this table, the data collected from the design using ODEs are related to the design of the sewage treatment plant system. The modeling was not based on ODE values. Instead, the ODE values were derived from figures provided by the World Health Organization (WHO) for WWTP. On the other hand, the observational data refer to the data gathered within a 15-day period from the actual performance of the sewage treatment plant. It is evident that the utilization of the SDE method has led to a reduction in model output uncertainties. Considering the utilization of ODEs in the design section to obtain a solitary outcome for each variable, average observational data and SDE over a span of 15 days were used to compare the results. The findings unequivocally demonstrate a good agreement between the results obtained through the SDE method and the observational data. The proposed method significantly improves the results (derived from Equations (16) and (17)) by considering the randomness of the problem, resulting in an enhancement of 11.47 and 10.11% for the BOD and *Q* models, respectively.

SDE (average) in 15 days . | ODE (designing value) . | Observation (average) in 15 days . | Variable . |
---|---|---|---|

22.834 | 20.00 | 23.0625 | BOD_{eff} (mg L^{−1}) |

142,707.5 | 129,600 | 142,660.9 | Q_{eff} (m^{3} day^{−1}) |

SDE (average) in 15 days . | ODE (designing value) . | Observation (average) in 15 days . | Variable . |
---|---|---|---|

22.834 | 20.00 | 23.0625 | BOD_{eff} (mg L^{−1}) |

142,707.5 | 129,600 | 142,660.9 | Q_{eff} (m^{3} day^{−1}) |

## CONCLUSIONS

Numerous physical and engineering phenomena exhibit a random nature, and neglecting the stochastic terms when modeling such processes with ODEs amplifies the uncertainty in the resulting outputs. Considering the assessment of uncertainty in modeling, SDEs offer a superior alternative to ODEs. This study aimed to derive the governing SDEs for the activated sludge process in the WWTP in Tabriz City, considering variations in the pollution parameters *Q*_{eff} and BOD_{eff} as biological operations. To address the fluctuations in the activated sludge process and numerically solve the integral from the SDE, the Euler–Maruyama method was utilized. The findings suggest that the biological aspect of WWTP can be explored through SDE. However, in order to avoid accumulating errors and account for the stochastic nature of the process, a 15-day interval was deemed optimal for determining the constant coefficients and variables. The results indicate that the use of SDE over extended periods of time with close-time data, due to fluctuations in constant coefficients, may not yield satisfactory outcomes and is not recommended for WWTP design. The SDE technique has demonstrated its reliability in monitoring the variations in *Q*_{eff} and BOD_{eff} within short time periods. This solution effectively managed the uncertainty associated with the dynamic conditions of death and birth, active microorganisms, (i.e., regulating the optimal food concentration for these bacteria in the secondary sedimentation unit). It is important to emphasize that the utilization of these microorganisms in the WWTP system leads to substantial overhead expenses for the beneficiaries. In contrast to black-box methods, the SDE method offers the ability to model stochastic processes with minimal data, making it suitable for monitoring the performance of the activated sludge sector with limited observational data and within a short timeframe. Future studies should explore alternative methods for solving SDEs and modeling activated sludge, comparing them with the results presented in this paper. Subsequently, the outcomes from these models can be combined to enhance the accuracy of the output.

## AUTHOR CONTRIBUTIONS

R. S. Z. rendered support in project administration, formal analysis, investigation, resources, data curation, and writing the original draft; V. N. conceptualized the whole article, supervised the article, and devised the methodology; M. S.-F. investigated the article, brought the resources, wrote the original draft, and wrote the review and edited the article; H. G. supervised the article and developed the methodology; C.-Q. K. supervised the article and devised the methodology.

## DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

*Simulation of Wastewater Treatment Plants Modeled by a System of Nonlinear Ordinary and Partial Differential Equations*

*Recursive Bayesian Inference on Stochastic Differential Equations*