ON MODELLING WATER QUALITY WITH STOCHASTIC DIFFERENTIAL EQUATIONS

Abstract

Based on biochemical kinetics, a stochastic model to characterize wastewater treatment plants and dynamics of river water quality under the influence of random fluctuations is proposed in this paper. This model describes the interaction between dissolved oxygen (DO) and biochemical oxygen demand (BOD), and is in the form of stochastic differential equations driven by multiplicative Gaussian noises. The stochastic persistence problem for the model of the system is analysed. Further, a numerical simulation of the stationary probability distributions of BOD and OD by approximations of the stochastic process solution is presented. These results have implications for the prediction and control of pollutants.

1 Introduction

Modelling and simulation of biochemical systems has become an area of intensive research in the past decades. The evolution of these system has been modelled by deterministic reaction rate equations. In general, nonlinear physics models have been intensively and numerously applied to uncover the biochemical complexities. Specifically, differential equation models for the quality of river water and a wastewater treatment plant have already been constructed. They simulate the effects and interactions between dissolved oxygen (DO) and biochemical oxygen demand (BOD), algae and others. Aspects of these models have been considered by different authors (see, for example, [3, 6, 12] and the references therein). Moreover, stochastic differential equation models for the evolution of DO and BOD were developed [1, 4, 8, 11, 13, 14]. Typically, such a system behaves stochastically rather than deterministically. The randomness in these models arise from a number of factors, such as the inherent variability and natural heterogeneity (for example, atmospheric conditions), measurement errors (for example, calibration), complexity (for example, complete mixing), lack of data (for example, of nonpoint sources) and others. Also, such models represent biological, chemical and physical processes and changes in such aquatic ecosystems. In this context, water quality dynamics is inevitably affected by environmental noises which can have significant effects on OD and BOD components [1, 4]. Overall, the objective of these studies is to determine the evolution of OD and BOD variables under the effects of such factors.

Motivated by the above research, in this paper, we consider a stochastic version of a water quality model under the influence of random fluctuations. It describes the interaction between BOD and DO, and is in the form of stochastic differential equations driven by multiplicative Gaussian noises. The aim of this study is to characterize the persistence of the long-time behaviour of the stochastic system. Such results enable us to predict the behaviour of water quality and the impacts of the random noises on the dynamics of the system. The paper is organized as follows. In Section 2, we introduce the stochastic differential equation model system. In Section 3, we carry out a computational analysis. This includes analytical and numerical results concerning the stochastic persistence of the long-time behaviour of the system. Section 4 concludes the paper.

2 The stochastic differential equation system

Consider the dynamics model describing water quality in terms of BOD and DO interaction with rates of reaction kinetics of nonlinear BOD degradation and DO depletion of Michaelis–Menten type [12], through a first-order system of differential equations, where the BOD and DO concentrations are denoted by $x_1(t)$ and $x_2(t)$ , respectively, and t is the time. Then, the BOD–DO interaction system is given by

(2.1) $$ \begin{align} \begin{aligned} \frac{d x_1(t)}{dt}= & -k_1\frac{x_2(t)}{k+x_2(t)}x_1(t)+q, \\ \frac{d x_2(t)}{dt}= &-k_2\frac{x_2(t)}{k+x_2(t)}x_1(t)+\alpha(x_{2_c}-x_2(t)), \end{aligned} \end{align} $$

where $k_1$ is the BOD decay rate, $k_2$ is the DO deoxygenation rate, k is the half-saturated oxygen demand concentration, $\alpha $ is the reaeration rate, $x_{2_c}$ is the oxygen saturation constant and q is the source term for BOD.

The overall dynamics of BOD and DO (that is, the model in equation (2.1)) under the influence of random fluctuations can be described by the stochastic differential equation system:

(2.2) $$ \begin{align} \begin{aligned} d x_1(t)&= \bigg(-k_1\frac{x_2(t)}{k+x_2(t)}x_1(t)+q\bigg)\,dt+\sigma_1 x_1 \, dW_1(t), \\ d x_2(t)&= \bigg(-k_2\frac{x_2(t)}{k+x_2(t)}x_1(t)+\alpha(x_{2_c}-x_2(t))\bigg)\,dt+\sigma_2 x_2 \, dW_2(t), \end{aligned} \end{align} $$

where $W_1(t)$ and $W_2(t)$ are independent Gaussian white noises defined on complete probability space $(\Omega , \mathcal {F},(\mathcal {F}_t)_{t\geq 0}, \mathbb {P})$ , and nonnegative constants $\sigma _1$ and $\sigma _2$ represent the noise intensities. For convenience, it is assumed that the variables $x_1$ and $x_2$ are dimensionless, that is, the equations describe changes in relative concentration size.

3 Analysis for stochastic persistence

To study the stochastic persistence and the impacts of random noise, we first discuss the dynamics of the system in equation (2.2) in the absence of noise by determining the stability of the solution. Therefore, we examine the equilibria and perform a phase plane analysis.

For the equilibrium points, the steady-state solution is considered. The steady-state equilibrium point $E_s=(x_{1_s}, x_{2_s})$ is given by

$$ \begin{align*} E_s=\bigg(\bigg(1+\frac{k}{x_{2_s}}\bigg)\frac{q}{k_1}, x_{2_c}-\frac{k_2 q}{\alpha k_1}\bigg),\, k_1,\, \alpha,\, x_{2_s}\neq 0. \end{align*} $$

This interior equilibrium point exists whenever $x_{2_c}>k_2q/\alpha k_1$ and $q\neq 0$ . To analyse the stability of this equilibrium point, we determine the characteristic equation of the associated Jacobian matrix of equation (2.1) given by

$$ \begin{align*} J(x_1, x_2)= \begin{pmatrix} -\dfrac{k_1 x_2}{k+x_2} & -\dfrac{k_1 k x_1 x_2}{(k+x_2)^2} \\[.3cm] -\dfrac{k_2 x_2}{k+x_2} & -\dfrac{k_2 k x_1 x_2}{(k+x_2)^2}-\alpha \end{pmatrix}. \end{align*} $$

Using the Routh–Hurwitz criterion [2], which guarantees that all the roots of the characteristic equation $\lambda ^2+a\lambda +b =0$ have negative real part if and only if $a>0$ and $b>0$ , and considering that the model parameters are positive, the coefficients of the characteristic equation of the Jacobian matrix $J(x_{1_s}, x_{2_s})$ are

$$ \begin{align*} \begin{split} a = & \frac{k_1 x_{2_s}}{k+x_{2_s}} +\frac{k_2 k x_{1_s} x_{2-s}}{(k+x_{2_s})^2}+\alpha>0 , \\ b = & \frac{\alpha k_1 x_{2_s}}{k+x_{2_s}}>0. \end{split} \end{align*} $$

Hence, $E_s$ is a stable equilibrium point. To verify that the model has a stable interior equilibrium $E_s$ , we take parameter values as in [6]: $k_1=0.15, k_2=0.5, q=0.2, x_{2_c}=9.5, \alpha =0.35, k=3$ . The resulting phase portrait is shown in Figure 1.

Figure 1 This figure shows the phase portrait of equation (2.1) for given parameter values.

Then, to study the stochastic persistence and the impact of random noise, the model system in equation (2.2) will be used.

3.1 Analytical results

In this subsection, we present some analytical results concerning the stochastic persistence for the model system in equation (2.2). For this, we follow the method of Gray et al. [5].

Theorem 1. For any given initial values $x_i(0)=x_{i_0}, i=1,2$ , the solution of the model in equation (2.2) obeys

(3.1) $$ \begin{align} \limsup_{t\rightarrow \infty} x_i^{-1}(t)\geq \xi_i \quad a.s., \end{align} $$

and

(3.2) $$ \begin{align} \liminf_{t\rightarrow \infty} x_i^{-1}(t)\leq \xi_i \quad a.s., \end{align} $$

where $\xi _1$ and $\xi _2$ are solutions of the following equations:

$$ \begin{align*} \begin{aligned} f_1(\xi_1)&= \frac{q}{\xi_1}-\frac{k_1 x_{2_s}}{k}-\frac{\sigma_1^2}{2} =0, \\ f_2(\xi_2)&=\frac{\alpha x_{2_c}}{\xi_2}- \frac{\alpha k+k_2 x_{1_s}}{k}- \frac{\sigma_2^2}{2}=0, \end{aligned} \end{align*} $$

respectively. (Here a.s. means almost surely.) That is, $x_i^{-1}(t)$ will rise to or above the level $\xi _i$ infinitely often with probability one.

Proof. We consider approximations for $dx_1(t)$ and $dx_2(t)$ in equation (2.2) and prove the assertion in equation (3.1). If it is not true, then there is a sufficiently small $\epsilon \in (0,1)$ such that

$$ \begin{align*} \mathbb{P}(\Omega_1)>\epsilon, \end{align*} $$

where $\Omega _1=\{\limsup _{t\rightarrow \infty } x_i^{-1}(t)\leq \xi _i-2\epsilon \}$ . Hence, for every $\omega \in \Omega _1$ , there is a ${T=T(\omega )>0}$ such that

(3.3) $$ \begin{align} x_i^{-1}(t,\omega) \leq \xi_i-\epsilon \quad \text{whenever } \, t\geq T(\omega). \end{align} $$

It follows from equation (3.3) that

(3.4) $$ \begin{align} f(x_i(t,\omega))\geq f(\xi_i-\epsilon) \quad \text{whenever } t\geq T(\omega). \end{align} $$

Moreover, by the large number theorem for martingales, there is a $\Omega _2\subset \Omega $ with $\mathbb {P}(\Omega _2)=1$ such that for every $\omega \in \Omega _2$ ,

(3.5) $$ \begin{align} \lim_{t\rightarrow \infty}\frac{1}{t}\int_0^t \sigma_i \, dW_i(s,\omega)=0. \end{align} $$

Now, fix any $\omega \in \Omega _1\cap \Omega _2$ . It then follows from the Ito formula [10] and equation (3.4) that, for $t\geq T(\omega )$ ,

$$ \begin{align*} \log (x_i(t,\omega))\geq \log(x_{i_0})+\int_0^{T(\omega)} f_i(x_i(s,\omega))\, ds +f_i (\xi_i-\epsilon)(t-T(\omega))+\int_0^t \sigma_i \,dW_i(s,\omega). \end{align*} $$

This yields

$$ \begin{align*} \liminf_{t\rightarrow \infty} \frac{1}{t} \log (x_i(t,\omega))\geq f_i(\xi_i-\epsilon)>0, \end{align*} $$

whence

$$ \begin{align*} \lim _{t\rightarrow \infty} x_i^{-1}(t,\omega)=\infty. \end{align*} $$

However, this contradicts equation (3.3). We therefore must have the desired assertion in equation (3.1).

Next we prove the assertion in equation (3.2). If it is not true, then there is a sufficiently small $\delta \in (0,1)$ such that

$$ \begin{align*} \mathbb{P}(\Omega_3)>\delta, \end{align*} $$

where $\Omega _3=\{\liminf _{t\rightarrow \infty } x_i^{-1}(t) \geq \xi _i+2\delta \}$ . Hence, for every $\omega \in \Omega _3$ , there is a ${\tau =\tau (\omega )>0}$ such that

(3.6) $$ \begin{align} x_i^{-1}(t,\omega)\geq \xi_i+\delta \quad \text{whenever } \, t\geq \tau(\omega). \end{align} $$

Now, fix any $\omega \in \Omega _3\cap \Omega _2$ . It then follows from the Ito formula that, for $t\geq \tau (\omega )$ ,

$$ \begin{align*} \log (x_i(t,\omega))\leq \log(x_{i_0})+\int_0^{\tau(\omega)} f_i(x_i(s,\omega))\, ds +f_i (\xi_i+\delta)(t-\tau(\omega))+\int_0^t \sigma_i \, dW_i(s,\omega). \end{align*} $$

This, together with equation (3.5), yields

$$ \begin{align*} \limsup_{t\rightarrow \infty} \frac{1}{t} \log (x_i(t,\omega))\leq f_i(\xi_i+\delta)<0, \end{align*} $$

whence

$$ \begin{align*} \lim _{t\rightarrow \infty} x_i^{-1}(t,\omega)=0. \end{align*} $$

However, this contradicts equation (3.6). We therefore obtain the desired assertion in equation (3.2). Hence, the proof is complete.

Moreover, we have that $x_i$ converges to a unique stationary distribution with the probability density given by

$$ \begin{align*} P_i(x_i)=\frac{b_i^{a_i}}{\Gamma(a_i)}x^{-(a_i+1)} e^{-b_i/x_i},\quad x_i>0,\quad i=1,2, \end{align*} $$

where

$$ \begin{align*} a_1&=\frac{2 \beta_1+\sigma_1^2}{\sigma_1^2},\quad b_1=\frac{2 q}{\sigma_1^2}, \quad \beta_1=\frac{k_1x_{2_s}}{k}, \\ a_2&=\frac{2 \beta_2+\sigma_2^2}{\sigma_2^2},\quad b_2=\frac{2 \alpha x_{2_c}}{\sigma_2^2},\quad \beta_2=\frac{k\alpha+k_2 x_{1_s}}{k}, \end{align*} $$

and $\Gamma (.)$ is the Gamma function, $\Gamma (x)=\int _0^\infty t^{x-1} e^{-t}\,dt$ , and we have

$$ \begin{align*} \lim_{t\rightarrow \infty} \frac{1}{t} \int_0^t x_i(s) \, ds\approx \zeta_i,\quad i=1,2, \quad \text{a.s.}, \end{align*} $$

respectively, with $\zeta _1=kq/k_1 x_{2_s}$ and $\zeta _2=k \alpha x_{2_c}/(k\alpha +k_2 x_{1_s})$ .

3.2 Computer simulations

Here, we present some computer simulations by applying the Milstein method [7, 9] to illustrate the behaviour of the model system in equation (2.2) and support the analytical results. The discrete equations are

(3.7) $$ \begin{align} \begin{aligned} x_{1_{j+1}} & =x_{1_j}+\bigg(-k_1\frac{x_{2_j}}{k+x_{2_j}}x_{1_j}+q \bigg)\Delta t \\ &\quad +\sigma_1 x_{1_j}\sqrt{\Delta t}W_{1_j}+\frac{\sigma_1^2}{2}x_{1_j}(W_{1_j}^2-1)\Delta t, \\ x_{2_{j+1}} & =x_{2_j}+\bigg(-k_2\frac{x_{2_j}}{k+x_{2_j}}x_{1_j}+\alpha (x_{2_c}-x_{2_j}) \bigg)\Delta t \\ &\quad +\sigma_2 x_{2_j} \sqrt{\Delta t}W_{2_j}+\frac{\sigma_2^2}{2} x_{2_j}(W_{2_j}^2-1)\Delta t, \end{aligned} \end{align} $$

where $\Delta t$ is time increment, and $W_{1_j}$ and $W_{2_j}\,(j=1,2,\ldots , n)$ are the independent Gaussian random variables.

We solve equation (3.7) with the time step $10^{-2}$ for computation of the stochastic solution. We consider the influence of noise for the same parameter values ${k_1=0.15, k_2=0.5, q=0.2, x_{2_c}=9.5, \alpha =0.35, k=3}$ , where the system in equation (2.1) possesses a stable equilibrium $E_s=(x_{1_s}, x_{2_s})=(1.89,7.52)$ . Resulting solutions for initial conditions $(x_1(0),x_2(0))=(12,6)$ are presented in Figures 2 and 3. In Figure 3, we show a plot of the sample paths of the solution $(x_1(t), x_2(t))$ for noise intensities $\sigma _1=\sigma _2= 0.1$ . Figure 3 also shows sample paths of $x_1(t)$ and $x_2(t)$ for large noise intensities $\sigma _1=\sigma _2= 0.4$ . We note that we obtain the same behaviour for different initial conditions. These figures clearly illustrate persistence and the effect of changing the noise intensities. As we can observe, for small noise intensities, the system shows small-amplitude noisy oscillations while for larger noise intensities, the system exhibits large-amplitude oscillations of complex form.

Figure 2 Sample paths of the system in equation (2.2) with small values of the noise intensities represented by red lines and its deterministic case represented by blue lines. Here, $\sigma _1=\sigma _2= 0.1$ . (Colour available online.)

Figure 3 Sample paths of the system in equation (2.2) with large values of the noise intensities and its deterministic case. Here, $\sigma _1=\sigma _2= 0.4$ . (Colour available online.)

The transition from small- to large-amplitude stochastic oscillations with increasing noise are accompanied by the significant changes of the form of the probability density. In Figures 4, 5, 6 and 7, we show the approximate probability density functions of the stationary distribution obtained by computer simulation of sample paths of $x_1(t)$ and $x_2(t)$ of equation (3.7) in the case of keeping the parameters the same but with different noise intensities $\sigma _1=\sigma _2= 0.1, 0.2, 0.3$ and $0.4$ , respectively. The simulations were run for $2$ million iterations with step size $10^{-2}$ . From these figures, one can see the changes of probability densities under increasing values of the noise intensities $\sigma _1$ and $\sigma _2$ . The probability densities reveal a maximum indicating the most probable concentration size and influence by the noise intensity. This maximum of the stationary distribution depends on the noise intensity. Clearly, the maximum is shifted for large noise intensities. These computer simulations of the solution also show that for large noise intensities, the distribution of the solution is skewed in the sense that the distribution of the solution becomes more symmetric for small noise intensities.

Figure 4 Plot of the probability density functions of the stationary probability distributions of the variables $x_1(t)$ and $x_2(t)$ for small noise intensities, corresponding to Figure 2, $\sigma _1=\sigma _2=0.1$ .

Figure 5 Plot of the probability density functions of the stationary probability distributions of the variables $x_1(t)$ and $x_2(t)$ for increasing the noise intensities $\sigma _1=\sigma _2=0.2$ .

Figure 6 Plot of the probability density functions of the stationary probability distributions of the variables $x_1(t)$ and $x_2(t)$ for increasing the noise intensities $\sigma _1=\sigma _2=0.3$ .

Figure 7 Plot of the probability density functions of the stationary probability distributions of the variables $x_1(t)$ and $x_2(t)$ for large noise intensities, corresponding to Figure 3, $\sigma _1=\sigma _2=0.4$ .

We note also that values of means of these probability distributions of $x_1(t)$ and $x_2(t)$ are $0.5266$ and $5.0379$ , respectively, which are useful in assessing the behaviour of the solution process.

4 Conclusion

One of the major concerns of many communities is to monitor pollution or water quality and its effects in various life processes in a stream. In general, to study the level of oxygen-related pollution, it has been agreed that the main assumption underlying the modelling approach is that the two variables, namely the concentration of BOD and that of DO, are sufficient to evaluate the quality of water. DO is a commonly used index of water quality since it reflects the general healthy state, or otherwise, of an aquatic environment. BOD, in company with measures such as suspended solids, ammonia and nitrate concentrations, characterizes the polluting load of the complex organic materials. Also, it is responsible for the removal of DO from the water, and hence the importance of identifying the dynamic relationships which govern the BOD–DO interaction. Therefore, many descriptions of the BOD–DO interaction have been proposed using mathematical modelling which are mostly based on a deterministic approach.

In this paper, we have studied a stochastic water quality model for persistence of the system under the effect of random fluctuations. Such a model describes the interaction between DO and BOD components, and is in the form of stochastic differential equations driven by multiplicative noises. We have analysed the persistence of the long-time behaviour of the system using techniques from dynamical systems and the theory of stochastic processes. Further, we have numerically simulated the stationary probability distributions of BOD and DO variables by approximating the solution process. These computer simulations assess the persistence of the long-time behaviour of the system. In particular, the stationary probability distribution has a maximum, indicating a very probable state. This maximum is, for instance, enhanced with decreasing noise intensities.

The implications of the study are briefly summarized as follows. First, the results obtained in this paper enable us to predict the behaviour of river water quality using DO and BOD variables, and the impacts of the random noises on the dynamics of the system. Second, this stochastic approach provides a practical method for a measurement of the risk to water quality through BOD loading capacity. Finally, these results are worthwhile for the prediction and control of pollutants.

In future work, such water quality models, which represent biological, chemical and physical processes, and describe nondegenerate diffusive systems, can be studied to consider more complex stochastic differential equations.