American Journal of Applied Mathematics
Volume 3, Issue 6, December 2015, Pages: 335-340

Modeling and Numerical Simulation of River Pollution Using Diffusion-Reaction Equation

Tsegaye Simon, Purnachandra Rao Koya*

School of Mathematical and Statistical Sciences, Hawassa University, Hawassa, Ethiopia

Email address:

(T. Simon)
(P. R. Koya)

To cite this article:

Tsegaye Simon, Purnachandra Rao Koya. Modeling and Numerical Simulation of River Pollution Using Diffusion-Reaction Equation. American Journal of Applied Mathematics. Vol. 3, No. 6, 2015, pp. 335-340. doi: 10.11648/j.ajam.20150306.24

Abstract: In the present study we have applied diffusion – reaction equation to describe the dynamics of river pollution and drawn numerical solution through simulation study. The diffusion-reaction equation is turn to be a partial differential equation since the independent variables are more than one that include spatial and temporal coordinates. The diffusion-reaction equation is widely applied to environmental studies in general and to river pollution studies in particular. River pollution models are special cases and are included in the broad area known as environmental studies. The diffusion – reaction equation is characterized by the reaction term. When the reaction term depends on the concentration of the contaminants then the original single diffusion-reaction equation will evolve to be a system of equations and this lead to analytical problems. The diffusion-reaction equations are difficult to solve analytically and hence we consider numerical solutions. For this purpose we first separate diffusion and reaction terms from the diffusion-reaction equation using splitting method and then apply numerical techniques such as Crank – Nicolson and Runge – Kutta of order four. These numerical methods are preferred because the systems of equations are solved accurately and efficiently. Detailed discussion of the results and their interpretations are included.

Keywords: River Pollution, Dissolved Oxygen, Biological Oxygen Demand, Diffusion-Reaction Equation, Splitting Method, Simulation Study, Crank – Nicolson Method, Runge – Kutta Method

1. Introduction

Advection – Diffusion – Reaction (ADR) equations are partial differential equations (PDEs) dependent on temporal and spatial coordinates. The ADR equations can be used to model mathematically a wide range of natural phenomenon and explain their dynamics with respect to time. The applications of the general advection-diffusion-reaction equations are wide and numerous. For instance ADR equations are used to solve pollutant transport models in scientific disciplines ranging from atmospheric studies through medical science to chemo taxis [15].

However, in the present study we have chosen the one-dimensional Streeter – Phelps equation which describing the river self-purification model as a concrete example. The Streeter – Phelps equation is applied to model the amount of dissolved oxygen (DO) in a stream after waste water is discharged into the stream. The Streeter – Phelps model describes the amount of pollutant downstream as the pollutant travels with the stream velocity in the direction of the stream flow. When a pollutant is added to or introduced into a water source then the dissolved oxygen of the water decreases to a minimum level and then gradually recovers and finally reaches a saturation level. Further, following the Lagrangian approach we reduce the advection-diffusion-reaction equation into diffusion – reaction equation [67].

The Diffusion-reaction systems are mathematical models. These models are used to explain how the concentrations of one or more substances are distributed in space and how these concentrations vary under the influence of the two processes viz., diffusion and reaction. Diffusion causes the pollutant substances to spread out in the river water and during the local chemical reactions the pollutant substances are transformed into each other [89].

The Diffusion-reaction systems are interesting on many levels, displaying phenomena such as pattern formation far from equilibrium, Turing structures, nonlinear waves such as solitons or spiral waves and spatial – temporal chaos. The efficient and accurate simulation of such systems, however, is difficult. This is because they couple a stiff diffusion term with a typically strongly nonlinear reaction term. When discretized this leads to large systems of strongly nonlinear, stiff ODEs. There has been much activity over the last few years in developing time stepping algorithms to deal with such problems in the numerical analysis community. This community demonstrated the most popular numerical algorithm currently being used to solve diffusion-reaction equations. That is a second order central difference scheme in space, coupled with an explicit forward Euler time stepping scheme. This seems an attractive method for two reasons. Firstly it is easy to implement, and secondly many people feel that there is an element of ’overkill’ when using a highly accurate high order method for a problem that does not require accurate solutions. It is unattractive, however, in that it is both inaccurate for studying spatial – temporal chaos, for example and inefficient. Moreover finite-difference methods can sometimes lead to fake solution [5]. In this paper, we used some of this method and Crank – Nicolson method for discretization of diffusion term and Runge – Kutta of order four methods to solve reaction term. In order to find the numerical approximation of the given problem, splitting methods is also used [1013].

2. Mathematical Model

Suppose that a polluted river contains  contaminants, with concentrations, where. Then a possible approach to model river purification to each of the chemical contaminants is given in (1) [1415].


The physical interpretations of various terms, variables and parameters used in (1) are as follows:represents advection,  represents diffusion, represent reaction,  measures distance along the direction of river flow,  represents river velocity, is cross – sectional area of the river, is concentration of contaminant, is net rate of addition of the suspension,  is diffusivity,  is the emission of the contaminant, is chemical reaction of contaminantwith other contaminants and is the time.

Since the chemical reaction depends on the concentrationsof the contaminants, equation (1) becomes coupled. This implies it is usually not possible to solve analytically. So to avoid this difficulty letbe the concentrations of the most important suspensions. Let be an appropriate measure of concentration of all other suspensions combined. Letbe the concentration of dissolved oxygen (DO) which is the most important variable in the purification of river and any substance that consume oxygen considered as pollutant since organisms underwater die without oxygen; and let be biochemical oxygen demand (BOD) which is the amount of oxygen the pollutants would need for their complete oxidization, per unit volume of river water.

For the model (1) we assumed that advection along the river be neglected or ignored that is; since we considered that the system is based on Lagrangian description. In the Lagrangian description the advection terms disappear, whereas they remain in the alternative Eulerian description. We assumed that pollution input has ceased and BOD can decay only by combining with oxygen or by flowing downstream; pollutants do not, for example, evaporate. Hence and also since oxygen is destroyed by chemical reactions. We have also taken and  In view of these assumptions the purification model (1) simplifies to the system of equations (2) and (3) below:



In the model described by (2) and (3),denotes dissolved oxygen diffusion coefficient and  denotes biochemical oxygen demand diffusion coefficient. To find suitable expressions for and the first and second order reaction rates and of the model equations (2) and (3), let us consider that oxygen diffuses into the river from the air immediately above the water and the air – water interface behaves like a membrane that is permeable to oxygen. Then the flux of oxygen into the river per unit area is given by. Hereis the effective depth of the imaginary membrane, is permeability to oxygen, and is the concentration of oxygen in the air immediately above the river surface. Up on multiplying the fluxwith the width of the river, we obtain the rate at which oxygen enters the river per time duration and is given in (4) below:


In (4), the parameteris a constant and has a dimension of a specific rate. The termrepresents the oxygen deficit and plays the role of are aeration coefficient. Let us now specify the reaction ratesandof the model described by (2) and (3). The chemical reaction that consumes dissolved oxygen of the river can be written symbolically as shown in (5) below:


Suppose that the rate at which the reactants, namely dissolved oxygen and biochemical oxygen demand, convert into product is proportional to their concentration. Moreover, by definition of biochemical oxygen demand, both the reactants must convert at the same rate. This conjunction results in a relation and is given in (6) below:


In (6), the parameteris a constant. The equation (6) also implies that and. Using equations (2), (3) and (4) and also introducing the function defined byin the equation (6) we obtain a pair of coupled non – linear system of partial differential equations given in (7) and (8) below:



Here in (7) and (8), the variable represents the spatial coordinate and the parameter represents the dissolved oxygen saturation level. The functionrepresents the deoxygenation and takes the values as given in (9) below:


In the equation (9), without loss of generality can be chosen as. The choice for in (9) guarantees that deoxygenation due to first and second order kinetics coincide when the river is fully saturated with dissolved oxygen, that is.

3. First Versus Second Order Kinetics

The vast majority of literature on environmental modeling [16 – 18] employs a first order kinetic model for the deoxygenation process. This has the sizable benefit of linearization of the problem and, if the diffusive effects are ignored, actually permitting an exact solution. However, we will now assume that if the river is even moderately polluted, then second order kinetics are in order.

4. Numerical Simulations

Let us solve the system of equations (7) and (8) with second order kinetics. As a test example, we make use of the system of equations given in (10) and (11) below:



Together with the system of equations (10) and (11), for the purpose of numerical simulation, let us restrict the independent variables to the regionsand. The boundary conditions are considered asand. The initial conditions are considered asand. Also, the parametric values are chosen asand.

4.1. Numerical Simulation for Zero – Diffusion

We begin by considering the case when the diffusion can be taken to be zero, i.e. . In this case the system of partial differential equations (10) and (11) becomes a system of ordinary differential equations given in (12) and (13) below. This corresponds to a well – mixed case. Thus, we are only studying the effect of self –purification, without considering the spatial distribution. Let the time interval of the integration is and the number of steps or subintervals in to which the time interval is divided be. Then the time durations required per one step is.



We now employ the method of Runge – Kutta of order – four on the system of equations given in (12) and (13). For this very purpose we consider the conditions and parametric values given just below the system of equations (10) and (11). The results of the simulation study are given in Figure 1.

Figure 1. Numerical simulation of the system of equations given in (12) and (13).

The carrying capacity of dissolved oxygen is considered to be one unit, i.e. . If at any time then the self-purification system of the water is not active. But, for any reason if at any time then the self-purification system of the water becomes active and helps the dissolved oxygen to boost up to reach its carrying capacity i.e. . As long as the dissolved oxygen does not reach its carrying capacity the self – purification system does not become inactive. But, when the dissolved oxygen reaches its carrying capacity, the self – purification system becomes inactive. Thus, the responsibility of the self – purification system of the water is to see always that the amount of dissolved oxygen be at or reach its carrying capacity or saturation level.

In Figure 1, the blue and red curves respectively represent the amounts of dissolved oxygen and the biochemical oxygen demandat any time in the water. At time the amount of dissolved oxygen and the amount of biochemical oxygen demand. Since the biochemical oxygen demandis positive, thetakes oxygen from. As a result the biochemical oxygen demand and dissolved oxygen availabilityare bothdecreasing. Now at this situation the self – purification capacity of the water becomes active. The water purification capacity helpsto do two things: (i) to reduceto zero and (ii) to increaseto reach its carrying capacity one unit. This scenario is pictorially described.

In other words, Figure 1 can be interpreted as follows: At, someone puts waste or pollutant in the water with biochemical oxygen demand concentration 6 times as high as dissolved oxygen concentration. The waste immediately reacts with the dissolved oxygen in the water causing the dissolved oxygen concentration to drop and also the biochemical oxygen demand concentration. After a while, the self-cleaning system of the water becomes active, so biochemical oxygen demand concentration goes down to zero and the dissolved oxygen concentration goes its normal value 1.

4.2. Numerical Simulation for Nonzero – Diffusion

We now consider that both the diffusion coefficients and  appear in the system of equations (10) and (11) are different from zero and see their effect through a simulation study. Note that in absence of the diffusion coefficients, i.e. , the system of equations (10) and (11) reduces to a system of ordinary differential equations. But in presence of the diffusion coefficients i.e. and the system of equations (10) and (11) remain to be a system of partial differential equations. In general the spatial axis is considered in the interval, i.e. . But, for the purpose of numerical simulation here we consider one unit length of spatial axis i.e. and hence.

We now discretize the spatial axis  into number of steps. Thus we have a step size. In general, the interval of thestep is given by. But in our present case it is. Applying the method of central difference in space, the system of equations (10) and (11) takes the form given by (14) and (15) below [19 – 22]:



In the system of equations (14) and (15), we have used the notationsand. Alsois a function of time. Let for Now, the system of equations (14) and (15) can be expressed in a matrix form shown in (16) and (17).



Here in (16) and (17), the symbol represents a tri – diagonal matrix of order  i.e.. Also the symbols, andare all vectors of dimension Further, represents an element by element product of the vectors and. To discretize the time coordinate, let us consider the integration on the time interval and divide the time interval into number of steps. Then the step size of time coordinate is given by. Up on splitting the diffusion equation from the reaction equation from the system of equations given by (16) and (17) we get a system of equations as shown in (18) and (19).



In short we can rewrite the system of two equations (18) and (19) into a single equation as, where  is an index, for diffusion (12) and reaction (13) equations respectively. Note that splitting diffusion from the reaction term has computational advantages since simultaneous coupling over space and the various chemical species is then avoided, and it also offers room for massively parallel computing [23 – 25]. Using step sizesolving the system of equations (18) and (19) by trapezoidal rule or Crank – Nicolson method and also solving the reaction equation (13) by Runge – Kutta of order four method, we have the simulation results with different diffusion coefficients as shown in Figure 2.

Figure 2. Numerical solution of (16) and (17) with and.

The Figure 2 represents concentrations of the dissolved oxygen and that of the biochemical oxygen demand. Figure 2 is obtained by the simulation study of the equations (16) and (17) with the conditions and parametric values given just below the system of equations (10) and (11). We have applied trapezoidal rule or also known as Crank – Nicolson method on the diffusion term and the method of Runge – Kutta of order four on the reaction term of the system of equations (16) and (17). In the case that the diffusion terms different from zero, we ignore all variation along the river and all variations in depth, and only look at a cross – section of a river of width one unit i.e. normalized. Again, at time t = 0, waste water is poured into the water, all over. But we assume that the water at the boundaries cleans. So there, the oxygen concentration is 1, and the BOD concentration is zero. But then, in addition to the self – cleaning effect from the reactive term, we get diffusion of the clean water from the boundaries, which makes the river cleans faster than if we did not have this effect.

Figure 3. Numerical solution of (16) and (17) with and.

Figure 4. Numerical solution of (16) and (17) with and.

In Figures 2 to 4 the profile for varying the diffusive term, we saw that when the rate of diffusion is high then the concentration of contaminant decreases faster and also the river cleans faster.

5. Conclusions

In this paper we have presented one example of advection-diffusion-reaction equation of the environmental or river water purification model and solved it by numerical methods. We considered the equation based on Lagrangian description. Since in Lagrangian description advection term disappears and the diffusion-reaction equation remains. A system of diffusion-reaction equations are coupled by the term of reaction. This system of equations (7) and (8) are decoupled when the first order kinetics are used and are coupled when second order kinetics are used. However, in this study we described numerical approximation techniques for the system of equations (7) and (8) by employing second order kinetics. The algorithm has been implemented in Matlab and generated the simulated graphics.


  1. David J. and Logan (2006). Applied Mathematics, John Wiley and Sons, Interscience.
  2. Tsegaye Simon (2013). Numerical Simulation of Diffusion – Reaction Equations: Application from River Pollution Model, Hawassa University, Hawassa, Ethiopia (Unpublished M. Sc. Thesis).
  3. Won Y., Wenwu C., Tae – Sang C. and John M. (2005). Applied Numerical Methods Using MATLAB.
  4. Aly – Khan K. (2003). Solving Reaction – Diffusion Equations 10 Times Faster.
  5. Sanderson A. R., Meyer M. D., Kirby R. M. and Johnson C. R. (2007). A Frame works for Exploring Numerical Solutions of Advection-Reaction-Diffusion Equations Using a GPU-Based Approach, Computer visual Sci. Vol. (10), PP.1-16.
  6. Holzbecker E. (2007), Environmental Modeling: Using MATLAB.
  7. Nas S. S., Bayram A., Nas E. and Bulut V. N. (2008). Effects of Some Water Quality Parameters on the Dissolved Oxygen Balance of Streams, Polish J. of Environ. Stud. Vol.17, PP. 531-538.
  8. Craster R.V. and Sassi R. (2006). Spectral Algorithms for Reaction – Diffusion Equations.
  9. Gerischa A. and Chaplain M.A.J. (2004). Robust Numerical Methods for Taxis – Diffusion – Reaction systems: Applications to Biomedical Problems. Mathematical and Computer Modeling 43 (2006), Pp. 4975.
  10. Hamdi A. (2006). Identification of Point Sources in Two Dimensional Advection – Diffusion – Reaction Equation: Application to Pollution Sources in a River. Stationary Case. Inverse Problems in Science and Engineering, 15, 8, 885–870.
  11. Scott A.S. (2012). A Local Radial Basis Function Method for Advection-Diffusion-Reaction Equations on Complexly Shaped Domains.
  12. Sportisse B. (2007). A Review of Current Issues in Air Pollution Modeling and Simulation, Computer Geo Sci. Vol. 11, Pp.159-181.
  13. Verwer J., Hundsdorfer G., Willem H. and Joke G. (2002). Numerical Time Integration for Air Pollution Models, Surv. Math. Ind. Vol. 10, Pp.107–174.
  14. Brain J. McCartin Sydney B. and Forrester Jr. (2002). A Fractional Step – Exponentially. Fitted Hopscotch Scheme for the Streeter – Phelps Equations of River Self – purification, Engineering Computations. Vol. 19(2), and Pp. 177–189.
  15. Mesterton – Gibbons M. (2007). A Concrete Approach to Mathematical Modelling, John Wiley and Sons.
  16. Kiely G. (1997). Environmental Engineering, McGraw-Hill.
  17. Mihelcic J.R. (1999). Fundamentals of Environmental Engineering, Wiley.
  18. Schnoor J. (1996). Environmental Modeling: Fate and Transport of Pollutants in Water, Air, and Soil, Wiley – Interscience.
  19. Evans G., Blackdge J. and Yardley. (2000). Numerical Methods for Partial Differential Equations. Springer – Verlag, London.
  20. KværnøA. (2009). Numerical Mathematics, Lecture Notes in TMA4215.
  21. Owren B. (2012). TMA4212 Numerical Solution of Partial Differential Equations with Finite Difference Methods.
  22. Shahraiyni Taheri H. and Ataie B. (2009). Comparison of Finite Difference Schemes for Water Flow in Unsaturated Soils, International Journal of Aerospace and Mechanical Engineering. Vol. 3(1), Pp. 1–5.
  23. Huangsdorfer W. (1996). Numerical Solution of Advection – Diffusion – Reaction Equations, Lecture Notes for a PhD Course, CWI Netherlands.
  24. Hundsdorfer W. and Verwer J.G. (2007). Numerical Solution of Time – Dependent Advection – Diffusion – Reaction Equations.
  25. Yazici Y. (2010) Operator Splitting Methods for Differential Equations.

Article Tools
Follow on us
Science Publishing Group
NEW YORK, NY 10018
Tel: (001)347-688-8931