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, [Reference Cunha, Coneglian and Poletti3, Reference Guaca and Poletti6, Reference Pimpunchat12] and the references therein). Moreover, stochastic differential equation models for the evolution of DO and BOD were developed [Reference Boano, Revelli and Ridolfi1, Reference Curi, Unny and Kay4, Reference Kallianpur and Xiong8, Reference Papadopoulos and Tiwari11, Reference Revelli and Ridolfi13, Reference Stijnen, Heemink and Ponnambalam14]. 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 [Reference Boano, Revelli and Ridolfi1, Reference Curi, Unny and Kay4]. 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 [Reference Pimpunchat12], 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
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:
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
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
Using the Routh–Hurwitz criterion [Reference Chang and Chen2], 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
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 [Reference Guaca and Poletti6]: $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.
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. [Reference Gray, Greenhalgh, Hu, Mao and Pan5].
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
and
where $\xi _1$ and $\xi _2$ are solutions of the following equations:
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
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
It follows from equation (3.3) that
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$ ,
Now, fix any $\omega \in \Omega _1\cap \Omega _2$ . It then follows from the Ito formula [Reference Øksendal10] and equation (3.4) that, for $t\geq T(\omega )$ ,
This yields
whence
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
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
Now, fix any $\omega \in \Omega _3\cap \Omega _2$ . It then follows from the Ito formula that, for $t\geq \tau (\omega )$ ,
This, together with equation (3.5), yields
whence
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
where
and $\Gamma (.)$ is the Gamma function, $\Gamma (x)=\int _0^\infty t^{x-1} e^{-t}\,dt$ , and we have
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 [Reference Higham7, Reference Kloeden and Platen9] to illustrate the behaviour of the model system in equation (2.2) and support the analytical results. The discrete equations are
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.
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.
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.
Acknowledgment
I would like to thank the anonymous referees whose suggestions improved this manuscript.