- Email: [email protected]
Analytical and numerical results of dissolved oxygen and biochemical oxygen demand in non-uniform open channel
Ecological Modelling 252 (2013) 11–22
Contents lists available at SciVerse ScienceDirect
Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel
Letter to the Editor a r t i c l e
i n f o
Keywords: BOD and DO Finite element method Dispersion Advection
a b s t r a c t A mathematical model for DO (dissolved oxygen) and BOD (biochemical oxygen demand) in a nonuniform open channel is presented in this work. This model can be used to study transport processes in one-dimensional flow. An accurate numerical scheme based on the finite element method with linear polynomial basis is developed for studying flows both upstream and downstream with respect to a point source. The accuracy of numerical solutions have been studied and compared with certain analytical solutions at steady state. Various test cases including only dispersion effect, DO sag curve, and the linear variation of cross-sectional area are investigated. We have also derived an analytical solution for the last case when the model is simplified under certain conditions. This model has been applied to real flow parameters of Tha Chin river in Thailand. Furthermore, the impacts of flow velocity, kinetic reaction order, and dispersion coefficients have been simulated to observe the dissolved oxygen and biochemical oxygen demand dynamics along the channel. Crown Copyright © 2012 Published by Elsevier B.V. All rights reserved.
1. Introduction Water pollution from human activities including social and economic development, is a crucial problem in environmental management. It directly affects human health and living quality. Also, pollutants from agricultural operations can play a dominant role in the degradation of surface and groundwater quality (Knight et al., 2000). The development of mathematical model enables us to gain better understanding of how to control and predict water quality in rivers or channels. Models generally require various water quality factors to be observed and used as inputs. For example, to model the pollutant distribution in a stream, the important indicators are the level of DO (dissolved oxygen) and BOD (biochemical oxygen demand) concentrations. Flow velocity and saturated oxygen concentration are also the key parameters of the model. The purpose of the model would be to predict the dissolved oxygen concentration and the biochemical oxygen demand concentration along the channel as well as the interaction of these two quantities. Recently, Pimpunchat et al. (2009) have presented a mathematical model of pollution in rivers. They have shown some analytical solutions for various situations of steady state flow. A model involving dissolved oxygen balance in rivers is proposed by Jha and Bhatia (2007), and this model has been improved in various ways. Historically, the famous mathematical model used to predict water quality in rivers was proposed by Streeter and Phelps (1925). They developed the model equations to study how much the declination of dissolved oxygen concentration in a river at certain given points downstream. At steady state flow, the analytical solution can be derived and represented in the form of DO sag curve solution. However, for some extended and complicated models, analytical solutions cannot be obtained. For such cases, efficient numerical methods are required to find approximate solutions of the model. Previously, various numerical methods have been presented for solving pollution models. A finite element method for solving onedimensional problems involving convection-dispersion process and sources of pollutant in river is presented by Onyejekwe and Toolsi (2001). Cockburn et al. (1989) proposed a discontinuous Galerkin method for solving the model equation in conservation law form. They constructed and analyzed a numerical method based on the class of TVB (total variation bounded) schemes. Tyagi et al. (2009) have presented an alternative numerical scheme for predicting the BOD variation in time at various successive distances from a wastewater outfall in a stream. The differential equations include advection, dispersion, and biochemical decay effects. The significance of stream self-purification is proposed by Wang et al. (1978). Water quality is indicated by the level of dissolved oxygen concentration, and the solubility of dissolved oxygen in water. Special emphasis is placed on the modeling of biological reactions and dissolved oxygen saturation concentration in fresh water. The BOD and DO sag models that describe water quality in rivers have been presented by Le et al. (2007). They provide an alternative approach for solving the first-order BOD equation in conjunction with water quality models. In their approach, the differential equation of dissolved oxygen concentration is solved analytically by the application of Laplace transform method. Very recently, multi-scale analysis has been proposed by Wu et al. (2011a) for deriving the analytical solution of multi-scale environmental dispersion in wetland flow for the determination of the dominant variation of contaminant concentration. The environmental dispersion of two-zone shallow wetland extended from single zone is analytically explored by Wu et al. (2011b) in terms of the longitudinal evolution of lateral mean concentration under environmental dispersion. In this work, we present a modified mathematical model for dissolved oxygen concentration interacting with biochemical oxygen demand in a non-uniform open channel flow. The model is in the same form as the equations presented in Pimpunchat et al. (2009). However, in this study, we have considered and developed an accurate numerical method for the case of arbitrary cross-sectional areas on both upstream and downstream flows. The numerical scheme presented is based on the finite element method with linear Lagrange basis functions. The accuracy of numerical solutions depends on the choice of basis function. The advantage of finite element method is that the 0304-3800/$ – see front matter. Crown Copyright © 2012 Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ecolmodel.2012.09.026
12
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Nomenclature P(x,t) X(x,t) Dp Dx l
v A q k1 k2 k ˛ S t n
concentration of BOD (kg m–3 ) concentration of dissolved oxygen (kg m–3 ) dispersion coefficient of pollutant (m2 day–1 ) dispersion coefficient of dissolved oxygen (m2 day–1 ) length of channel (m) average of water velocity (m day–1 ) cross-sectional area (m2 ) rate of pollutant addition (kg m–1 day–1 ) degradation rate coefficient for pollutant (day–1 ) de-aeration rate coefficient for dissolved oxygen (day–1 ) half-saturated oxygen demand concentration for pollutant decay (kg m–3 ) mass transfer of oxygen from air to water (m2 day–1 ) saturated oxygen concentration (kg m–3 ) time (days) number of elements in the finite element method
scheme can be extended directly to higher order by changing the order of basis approximation. In contrast, in the finite difference method derivatives are approximated in terms of difference equations locally and usually suffer from poor approximation of nonlinear term. Using the finite element method, the concentration of dissolved oxygen interacting with biochemical oxygen demand can be approximated at specified positions and time. We have performed various test cases to validate the accuracy of our numerical scheme. Test cases include only dispersion effect, DO sag curve, linear variation of cross-sectional area, and the application to Tha Chin river in Thailand. We have also given an analytical solution for a simplified model. Furthermore, the effects of dispersion, advection, and reaction kinetics are studied by including these terms in the model, and then solving numerically. The mathematical model of dissolved oxygen concentration interacting with biochemical oxygen demand concentration in an open channel is presented in Section 2. The finite element method and time marching scheme are presented in Section 3. Test cases and numerical results are shown in Sections 4 and 5 respectively. Finally, we have made some conclusions in Section 6.
2. Mathematical model We consider a one-dimensional flow along a channel. The spatial variable x is used to describe the distance along the channel. Environmental quantities, such as BOD and DO concentrations, are allowed to vary along the length of the channel. The mathematical model for the interaction between the BOD and DO concentrations along the channel includes the effects of advection, dispersion, and reaction. The model equations can be derived from the conservation of mass for an arbitrary cross-sectional area of the channel. The governing equations can be expressed as follows. ∂2 (AP) X ∂ (vAP) ∂ (AP) = Dp − K1 AP + qH(x), − X +k ∂t ∂x ∂x2
(1)
∂2 (AX) X ∂ (vAX) ∂ (AX) = Dx − K2 AP + ˛(S − X). − X +k ∂t ∂x ∂x2
(2)
This model includes dispersion and advection effects for unsteady processes in arbitrary cross-sectional area which is a function of distance. A reaction term is included in the form of Michaelis–Menten kinetics. The model allows adding pollution which is represented by the Heaviside function H. The reaeration process of dissolved oxygen until reaching the saturated concentration S is represented by the last term in Eq. (2). In this work, we have generalized the model from Pimpunchat et al. (2009) by allowing the cross-sectional area A to be a function of x. In some specific case of flows, analytical solutions are possible to derive. But, for more complex flows, analytical solution cannot be obtained. Thus, we apply the finite element method to solve numerically the model equations. The accuracy of numerical solutions is tested by comparing with some existing analytical solutions for some flow cases.
3. Finite element method The main steps of the finite element formulation for solving Eqs. (1) and (2) are presented as follows.
3.1. Spatial discretization We begin by dividing the computational domain (–l ≤ x ≤l) in n elements, element length is L = 2l/n. The points at x = –l and x = l refer to the end points of upstream and downstream respectively.
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
13
3.2. Approximate function The BOD concentration for each element is approximated by P (x) =
2
i Pi = 1 P1 + 2 P2 = [1 2 ]
i=1
P1
= [] [P] .
P2
(3)
Similarly, the dissolved oxygen concentration is approximated by
X (x) = [1 2 ]
X1
= [] [X] ,
X2
(4)
where 1 = 1 − x/L and 2 = x/L are the linear Lagrange basis functions for each element, [P], and [X] are the values of the unknown functions at each node of the element. 3.3. Weighted residual method The system (1) and (2) is approximated by the standard Galerkin finite element method based on weighted residuals. The system can be written as
L W
∂ (AP) ∂2 (AP) X ∂ (vAP) − Dp + K1 AP − qH(x) + X +k ∂t ∂x ∂x2
dx = 0,
(5)
0
L W
∂ (vAX) ∂ (AX) ∂2 (AX) X − Dx + K2 AP − ˛(S − X) + X +k ∂t ∂x ∂x2
dx = 0,
(6)
0
where W is a weighting function. We choose the weighting function to be the same as the basis function. Thus, the finite element formulation for Eqs. (5) and (6) becomes
L
∂ (AP) i dx − ∂t
L
∂2 (AP) i Dp dx + ∂x2
L
∂ (vAP) dx + i ∂x
L
0
0
0
0
L
L
∂ (AX) i dx − ∂t
∂2 (AX) i Dx dx + ∂x2
0
0
L
X APdx − i K1 X +k
0
i qH(x)dx = 0,
(7)
i ˛(S − X)dx = 0,
(8)
0
L
∂ (vAX) i dx + ∂x
L
X i K2 APdx − X +k
0
L 0
where i refers to the node number in an element (i = 1,2). 3.4. Finite element formulation For the special case when the half-saturated oxygen demand for pollutant decay, k, is negligible (k≈0), Eqs. (7) and (8) simplify to
⎡L ⎤ ⎡ L ⎤ ⎡ L ⎤ dj di dj ⎣ A(x)i j dx⎦ P˙ j + ⎣Dp A(x) dx⎦ Pj + ⎣v A(x)i dx⎦ Pj ⎡
0
+ ⎣K1
dx dx
L
0
⎤
A(x)i j dx⎦ Pj −
0
dx
0
L
,
(9)
,
(10)
qH(x)i dx = 0 0
⎡L ⎤ ⎡ L ⎤ ⎡ L ⎤ dj di dj ⎣ A(x)i j dx⎦ X˙ j + ⎣Dx A(x) dx⎦ Xj + ⎣v A(x)i dx⎦ Xj 0
⎡ + ⎣K2
L
⎤0
dx dx
A(x)i j dx⎦ Pj − ˛S
0
where P˙ j ≈ (Pjn+1 − Pjn )/t, Pj = X˙ j ≈ (Xjn+1 − Xjn )/t, Xj =
L
⎡
L
i dx + ⎣˛
0
Pjn+1 + (1 − )Pjn
0
⎤
dx
i j dx⎦ Xj = 0
0
and
(11)
Xjn+1 + (1 − )Xjn .
(12)
14
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Here Pj and Xj denote the variable’s values at time t = tn and at node xj , t is the time step, and is usually specified in the range of 0 ≤ ≤ 1 and is used to weight the approximate solutions in time. The most commonly used value of is 0 (gives the explicit forward Euler scheme), 0.5 (gives a second-order central implicit method, or Crank-Nicolson-Galerkin), or 1 (leads to the fully implicit method). In an element, we allow the cross-sectional area A to vary as a function of x representing arbitrary cross-sectional area. Thus we approximate A by A = [][A] where [A]T = [A1 A2 ] represents the cross-sectional area at two nodes of the element. The resulting Eqs. (9) and (10) can then be written in matrix forms as
[T ] + t
Ap + [B] + K1 [T ]
Pjn+1 = [T ] + t − 1
Ap + [B] + K1 [T ]
Pjn + t [G] ,
[T ] + t ([Ax ] + [B] + [F2]) Xjn+1 = [T ] + t − 1 ([Ax ] + [B] + [F2]) Xjn − t [[E] − [F1]] ,
where
L
L A(x)i j dx = 12
[T ] =
3A1 + A2
A1 + A2
A1 + A2
A1 + 3A2
0
(13) (14)
,
(15)
and dispersion matrices shown by
L [AP ] == Dp
Dp d dj dx = A(x) i 2L dx dx
0
and,
Dx [Ax ] = 2L
− (A1 + A2 )
− (A1 + A2 )
A1 + A2
A1 + A2
− (A1 + A2 )
− (A1 + A2 )
A1 + A2
The matrix of BOD advection is
L [B] = v
A1 + A2
A(x)i
dj dx
dx =
,
(16)
.
(17)
v
− (2A1 + A2 )
2A1 + A2
6
− (A1 + 2A2 )
A1 + 2A2
0
(18)
The matrix of adding pollution is given by [G] = qH(x)
L 2
1 1
.
(19)
Similarly, matrices [F1], [F2], and [E] are
L i dx = ˛S
[F1] = ˛S 0
1
i j dx =
˛L 3
⎢ ⎣
0
L A(x)i dx = K2
[E] = K2
,
1
⎡
L [F2] = ˛
L 2
0
(20)
1
1 2
1 2
1
L
⎤ ⎥ ⎦,
(21)
xA2 xA1 + A1 − L L
⎡ x⎤
1− L⎦ ⎣ dx.
(22)
x L
0
Finally, after assembling all elements together, we obtain the full form of the matrix system,
[T ]sys + t
sys
+ [B]sys + K1 [T ]sys
sys
+ [B]sys + K1 [T ]sys
Ap
= [T ]sys + t − 1
Ap
= [T ]sys + t − 1
Pjn+1
[T ]sys + t [Ax ]sys + [B]sys + [F2]sys
[Ax ]sys + [B]sys + [F2]sys
(23)
Pjn + t[G]sys ,
Xjn+1
Xjn + t[F1]sys − t[E]sys Pjn+1 − t(1 − )[E]sys Pjn .
(24)
Subscript sys means the ensemble matrices. These two systems can be solved iteratively by the Gauss-Seidel method. For all our computations, the tolerance is set to 10–9 in order to ensure completely converged solutions. The system (23) is solved iteratively before setting the values into the RHS of system (24) to obtain the approximation of Xjn+1 and Pjn+1 at the same time level. When we consider the interaction between BOD and DO, k cannot be neglected. The finite element systems (23) and (24) can be applied, but must be modified slightly. In this case, the Trapezoidal rule is employed for approximating the integral values in matrices [T]sys and [E]sys at time tn .
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
15
3.5. Case of variable dispersion For arbitrary cross-sectional-area as a function of distance, the dispersion coefficient can also be regarded as a function of distance. In this case, the governing equations become, ∂ (AP) ∂ = ∂t ∂x ∂ ∂ (AX) = ∂t ∂x
∂P Dp A ∂x
Dx A
∂X ∂x
−
−
X ∂ (vAP) − K1 AP + qH(x), X +k ∂x
∂ (vAX) X − K2 AP + ˛(S − X). X +k ∂x
Here Dp and Dx are functions of x. To obtain the finite element system, we approximate the dispersion coefficients by
Dp = [1 2 ]
Dp1
and
Dp2
Dx = [1 2 ]
Dx1
,
Dx2
where Dpi and Dxi for i = 1,2 are the dispersion coefficients of pollutant and DO at nodes of the element. The dispersion matrices [Ap ] and [Ax ] in Eqs. (13) and (14) are reformulated as
L [Ap ] =
Dp (x)A(x)
di dj dx, dx dx
Dx (x)A(x)
di dj dx. dx dx
0
L [Ax ] = 0
After integrating, we obtain the modified dispersion matrices for the case of variable dispersion as
Ap
(2Dp1 A1 + Dp2 A2 + Dp2 A1 + 2Dp2 A2 ) = 6L
[Ax ] =
(2Dx1 A1 + Dx2 A2 + Dx2 A1 + 2Dx2 A2 ) 6L
1
−1
−1
1
1
−1
−1
1
,
.
It can be seen that, if Dp1 = Dp2 and Dx1 = Dx2 , it means that Dx and Dp are constants. These resulting matrices are the same as those shown previously in Eqs. (16) and (17). So, for the case of variable dispersion coefficients and variable cross-sectional area, these modified matrices can be applied in place of Eqs. (16) and (17) for flow simulations. 3.6. Convergence of iterative scheme We solve the linear systems (23) and (24) by iterative scheme. So the convergence of our numerical solution depends on the properties of matrix coefficient. If we define, [PO] = [T ]sys + t
Ap
sys
+ [B]sys + K1 [T ]sys ,
and
[DO] = [T ]sys + t [Ax ]sys + [B]sys + [F2]sys ,
the convergence of the linear system (23) and (24) by the Guass–Seidel method depends on the invertibility of matrices [PO] and [DO]. The spectral radius, , which is defined as the largest magnitude of the eigenvalues of the matrix, can be used to analyze the convergence of the iterative scheme. To obtain a converged solution, we require the spectral radius of matrix [PO] denoted by PO , and matrix [DO] denoted by DO less than unity. Deriving explicit forms of stability criteria from these two matrices is difficult and needs to be analyzed in detail further. However, it can be seen that, if the element size is fixed, convergence depends directly on the values of t and , which are the time step and the implicitness of iterative scheme employed in the approximations. 4. Test cases In this section, we present numerical examples for investigating the efficiency of the finite element method. Three test cases are presented: only dispersion effect, the BOD–DO sag-curve behavior, and the linearly varied cross-sectional area. New analytical solutions are derived in this work for checking the accuracy of the numerical method. 4.1. Dispersion effect We consider the case when the dominant effect comes from dispersion that results from different levels of concentration in the longitudinal direction of flow. For simplicity, dispersion coefficient is taken as a constant, and so the exact solution is easily derived. Note that, for more realistic problems, dispersion coefficient should be a function of position, but exact solutions are not available for most problems. / 0 and However, the numerical scheme presented in section 3.5 can be applied for such cases. When k is assumed to be negligible (k≈0), q =
16
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
A is constant, the analytical solutions of BOD and DO concentrations at steady-state on the upstream and downstream flows are presented in Pimpunchat et al. (2009). The solutions are expressed as follows.
⎧ q ı+ˇ ⎪ ı+ˇ)x ( ⎪ 1 − e ,x ≥ 0 ⎪ ⎨ K1 A 2ˇ P(x) = , ⎪ ı+ˇ q ⎪ ı+ˇ)x ( ⎪ e ,x<0 ⎩ K1 A
(25)
2ˇ
v2 +4Dp K1
where ı = 2Dv , ˇ =
2Dp
p
, and
⎧ K2 q ı+ˇ + K2 q ı−ˇ ı + ˇ (ı+ˇ)x ⎪ (+ )x ⎪ + − S − + e − xe ,x ≥ 0 ⎪ ⎨ K1 ˛ K1 2 ˛ 4ˇ A∗ 4ˇ B∗ 2ˇA∗ X(x) = , ⎪ K2 q ı−ˇ ı − ˇ (ı+ˇ)x ı+ˇ − ⎪ (+ )x ⎪ + + e − xe , x < 0 − S ⎩ K1
4ˇ A∗
2 ˛
v2 + 4˛Dx
4ˇ B∗
(26)
2ˇB∗
A , A∗ = 2ADx ı − ˇ − vA , and B∗ = 2ADx ı + ˇ − vA. where = 2Dv x , = 2Dx To compare with the analytical solutions at steady state, we set the initial and boundary conditions in our simulations as
q P(0, t) = K1 A
1−
ı+ˇ 2ˇ
P(x, 0) = 0, X(x, 0) = S,
,
X(0, t) = S −
K2 q K2 q + K1 ˛ K1
+ ı−ˇ ı+ˇ + − 2 ˛ 4ˇ A∗ 4ˇ B∗
,
∂X(l, t) ∂P(−l, t) ∂X(−l, t) ∂P(l, t) = 0, = 0, = 0, and = 0 for t > 0 ∂x ∂x ∂x ∂x
4.2. BOD–DO sag curve When the half-saturated oxygen demand concentration for pollutant decay is negligible (k≈0), and there is no pollution added along the channel, except only the pollutant discharge point at x = 0, the analytical solutions for BOD and DO concentrations on the upstream and downstream flows can be obtained. The solutions at steady state are expressed as follows.
P(x) =
P0 e 2 x , x ≥ 0
(27)
P0 e 1 x , x < 0,
⎧⎛ ⎞ ⎪ ⎪ K K2 P0 P ⎪ 2 0
2 x ⎪ ⎝X0 − S − ⎠ e 4 x + + S, x ≥ 0 ⎪ ˛ ˛e ⎪ 2 2 ⎨ Dx 2 − v 2 − Dx 2 − v 2 − A A ⎛ ⎞ X(x) = ⎪ ⎪ P K K2 P0 ⎪ 2 0
1 x ⎪ ⎝X0 − S − ⎠ e 3 x + + S, x < 0, ⎪ ˛ ˛e ⎪ 2 2 ⎩ Dx 1 − v 1 − Dx 1 − v 1 − A
where
1 =
4 =
v+
v2
+ 4K1 Dp
2Dp
v−
v2 + 4 ˛A Dx 2Dx
,
2 =
v−
A
" v2
+ 4K1 Dp
2Dp
,
3 =
v+
˛ A
v2 + 4 Dx
, and
2Dx
.
To compare numerical results with the analytical solutions at steady state, we set initial and boundary conditions as P(x, 0) = 0, X(x, 0) = S, P0 = P(0, t) = 5, X0 = X(0, t) = S, ∂X(l, t) ∂P(−l, t) ∂X(−l, t) ∂P(l, t) = 0, = 0, = 0, and = 0 for t > 0. ∂x ∂x ∂x ∂x
(28)
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
17
4.3. Linear variation of cross-sectional area In this case, we assume that pollutant is being added along the non-uniform channel at constant rate of q = 0.2 kg m–1 day–1 . The mass transfer of oxygen from air to water ˛ is negligible. The reaction constant is zero (k = 0). To show the ability of our numerical scheme for solving the model in the case of non-uniform cross-sectional area, the function of A is simply assumed to be
# # ⎧ #b − a# ⎪ ⎪ ⎨ a − 200 x, x ∈ [0, 200] A(x) = # # ⎪ #b − a# ⎪ ⎩ a+
200
x, x ∈ [−200, 0) ,
where a and b are the values of cross-sectional area at x = 0 and x = 200 respectively. After solving the governing Eqs. (1) and (2), the analytical solutions can be obtained as follows.
⎧ q q A0 P0 − e 2 x + ⎪ ⎪ K K ⎪ 1 1 ⎪ ,x≥0 ⎨ A P (x) =
⎪ q q ⎪ ⎪ A0 P0 − e 1 x + ⎪ ⎩ K1 K1
X (x) =
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨
A
, x < 0,
q K1 Dx 22 − v 2
K2 A0 P0 − A0 X0 −
(29)
q K1 Dx 22 − v 2
K2 A0 P0 − +
e 2 x
A
,x ≥ 0
⎛ q ⎞ q ⎪ ⎪ K2 A0 P0 − K2 A0 P0 − ⎪ ⎪ K1 ⎠ 5 x K1 ⎪ e + e 1 x ⎪ ⎝A0 X0 − ⎪ 2 − v
2 − v
⎪ D
D
x x 2 2 ⎪ 2 2 ⎪ ⎩ , x < 0,
(30)
A where 5 = Dvx , A0 P0 = A(0,t)P(0,t), and A0 X0 = A(0,t)X(0,t) for t > 0. Initial and boundary conditions in our simulations are given by P(x, 0) = 1, X(x, 0) = S, P0 = P(0, t) = 1, X0 = X(0, t) = S, ∂X(l, t) ∂P(−l, t) ∂X(−l, t) ∂P(l, t) = 0, = 0, = 0, and = 0 for t > 0. ∂x ∂x ∂x ∂x 5. Numerical results In all of our simulations, we set the degradation rate coefficient as K1 = 0.5 day–1 , the de-aeration rate coefficient as K2 = 0.05 day–1 , the saturated oxygen concentration as S = 2 kg m–3 , the water velocity as v = 2 m day–1 , the dispersion coefficient of pollution or BOD as Dp = 5 m2 day–1 , the dispersion coefficient of dissolved oxygen as Dx = 5 m2 day–1 , and the final time of simulations as t = 200 days. This duration is very long, and is used to investigate the convergence of steady solution. By fixing these model parameters, we then simulate the test cases by varying pollutant added along the channel (q), and the mass transfer of oxygen from air to water (␣). 5.1. Dispersion effect We set model parameters to v = 2 m day–1 , q = 0.2 kg m–1 day–1 (for x≥0), q = 0 kg m–1 day–1 (for x<0), Dp = 5 m2 day–1 , Dx = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 . We consider the uniform cross-sectional area in this case. The results of BOD concentration on the upstream and downstream flows at various time steps are shown in Fig. 1. The initial BOD concentration is zero along the channel. The numerical solution of BOD concentration increases dramatically downstream and approaches the analytical solution described in Eq. (25) when k≈50 days. There is no pollution added to the upstream flow, but there is a small dispersion area near the point of x = 0. This shows dispersion effect from the pollution added on downstream. The long-time numerical solutions are in good agreement with the analytical solution at steady state. This shows the accuracy of the present scheme. Moreover, we can apply this numerical method to predict the dynamics of BOD concentration before approaching the steady state in which analytical solutions are not available. The numerical results for the dissolved oxygen concentration are shown in Fig. 2. The initial concentration is given by X(x,0) = S = 2. As time increases, dissolved oxygen concentration decreases exponentially on the downstream channel. Its value, at t≈70days, approaches the analytical solution described in Eq. (26) at steady state. The dissolved oxygen concentration decreases due to the effect of constant pollution adding along the downstream channel. It has no effects to the dissolved oxygen concentration on the upstream flow, except at small area near x = 0. 5.2. BOD–DO sag curve Initially, BOD concentration is zero along the uniform channel, except only at the discharge point of x = 0. At this point, we assume that the BOD concentration is released continuously described by P(0,t) = 5 t > 0. It corresponds to a continuous adding of a point source. There
18
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Fig. 1. BOD concentration versus length when v = 2 m day–1 , q = 0.2 (for x ≥ 0), q = 0 kg m–1 day–1 (for x < 0), Dp = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 .
is no pollutant added on the upstream and downstream flows. The dissolved oxygen concentration on the stream can be increased by the re-aeration process where we set the re-aeration constant as ˛ = 0.01 m2 day–1 . The numerical results of BOD concentration and the analytical solution at steady state are shown in Fig. 3. The unsteady numerical solutions show the dynamics of BOD concentration at various time steps. The BOD concentration gradually increases from the point source and reaches the steady value when t≈9 days. Influent zone downstream (x≈20) is larger than upstream (x≈2) due to the flow direction effect, since we set the flow from left to right. The numerical results of dissolved oxygen concentration and the analytical solution at steady state are shown in Fig. 4. The initial DO concentration is set as X(x,0) = S = 2 which is a constant value along the channel. Numerical simulations show that dissolved oxygen concentration decreases as time increases and approaches the steady value for very large time. DO concentration decreases rapidly in the downstream flow due to the point source effect and then increases again on the far downstream due to re-aeration process. The profiles of DO concentration on the downstream are called the DO sag curve. This simulation shows the effect of pollutant addition and self-recovery area on the downstream. 5.3. Linear variation of cross-sectional area For simplicity in order to derive analytical solution, we assume that the cross-sectional area is varied as a linear function ranging from 10 = m2 to 200 = m2 (for x = 0 to x = l) and from 200 = m2 to 10 = m2 (for x = −l to x = 0). The initial BOD and DO concentrations are set as 1 kg m–3 and 1.6 kg m–3 respectively. Other parameters are given by v = 2 m day–1 , Dp = 5 m2 day–1 , and Dx = 5 m2 day–1 . Pollutant is added on both upstream and downstream by specifying q = 0.2 kg m–1 day–1 , except only at the point of x = 0, where q is zero. There is no re-aeration process in this case, so we set ˛ = 0.
Fig. 2. DO concentration versus length when v = 2 m day–1 , q = 0.2 kg m–1 day–1 (for x ≥ 0), q = 0 kg m–1 day–1 (for x < 0), Dx = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 .
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Fig. 3. BOD concentration versus length when v = 2 m day–1 , Dp = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 .
Fig. 4. DO concentration versus length when v = 2 m day–1 , Dx = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 .
Fig. 5. BOD concentration versus length when v = 2 m day–1 , q = 0.2 kg m–1 day–1 , Dp = 5 m2 day–1 , and ˛ = 0 m2 day–1 .
19
20
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Fig. 6. DO concentration versus length when v = 2 m day–1 , q = 0.2 kg m–1 day–1 , Dx = 2 m2 day–1 , and ˛ = 0 m2 day–1 .
The numerical results of BOD concentration at various time steps are shown in Fig. 5. Due to the constant pollutant loading along the non-uniform channel, the BOD concentration decreases rapidly on the upstream flow near x = 0. The BOD concentrations further decrease and approach zero at the far upstream. Similar results appear downstream. These results reflect the increasing cross-sectional area effects in those positions. The numerical results of dissolved oxygen concentration and the exact solution at steady state from Eq. (30) are shown in Fig. 6. The constant DO is specified continuously at the point of x = 0 meaning the addition of DO concentration. For the same time, constant BOD is also specified at this point. Numerical simulations show that dissolved oxygen concentration decreases due to the constant point source addition both upstream and downstream along the channel. But, in the early downstream flow, the DO concentration does not decrease rapidly like the DO sag curve case. The rate of decrease is smaller due to the cross-sectional area effect that is represented as an increasing function on the left and right of the stream flow. 5.4. Simulation of Tha Chin river In this section, we apply the numerical scheme to simulate the change of DO and BOD concentrations in the Tha Chin river located in Thailand. The variables and flow parameters of this river are reported in Pimpunchat et al. (2009) and PCD (2000). In SI units, the length of river is 325 km, the water velocity is 43,200 m day–1 , the cross-sectional area is 2100 m2 , the rate of pollutant addition along the river is 0.06 kg m–1 day–1 , K1 = 8.27 day–1 , K2 = 44.1 day–1 , ˛ = 16.5 m2 day–1 , and S = 0.01 kg m–3 . We simulate the case of zero dispersion. Graphs of P and DO concentrations along the river are shown in Figs. 7 and 8 respectively. Distance in these figures is shown in units of kilometer, and the unit of concentration is kilogram per hour. We set DO concentration at initial time to 10 mg l–1 . If we add pollutant with constant value along the river, the DO concentration decreases dramatically along the river and approaches the steady solution. For approximately 5 days,
Numerical value of pollution concentration (P) and exact solution at steady state 3500 t = 20 hrs 3000
pollution (P)
2500
2000
1500
1000
500 t = 10 hrs 0
0
50
100
150
200
250
300
350
length of river (x) Fig. 7. BOD concentration versus length of Tha Chin river, dashed line is simulated results and solid line is the exact solution at steady state.
Letter to the Editor / Ecological Modelling 252 (2013) 11–22 6
10
x 10
21
Numerical value of dissolved oxygen (X) and steady solution t = 10 hrs
dissolved oxygen (X)
9.8
9.6
9.4 t = 120 hrs 9.2
9
8.8
0
50
100
150
200
250
300
350
length of river (x) Fig. 8. DO concentration versus length of Tha Chin river, dashed line is simulated results and solid line is the exact solution at steady state.
DO concentration drops from 10 to 9.3 mg l–1 at a distance of 200 km. It can be seen further that DO concentration remains decreased along the downstream river when pollutant added downstream. This shows that our simulator can be used to predict and give some viewpoints of pollution control in rivers. 6. Conclusions In this paper, we present a mathematical model describing the interaction between dissolved oxygen and biochemical oxygen demand in a non-uniform open channel. The model includes both upstream and downstream flows where the separating point can be regarded as a point-source location. Transport processes are assumed to be unsteady in one-dimensional space while flow velocity is constant. This model can also be used to study the dynamics of BOD and DO concentrations when pollution is added along the channel. DO is increased by the re-aeration process until approaching saturated concentration. Thus the solutions from this model provide some viewpoints for ecological management and pollution control in channels. To solve the model, we have developed an accurate scheme based on the finite element method. Numerical solutions in each element are approximated by applying the linear Lagrange basis functions. This then results in a coupled linear system of BOD and DO concentrations at arbitrary partitions of space and time. The accuracy of our numerical scheme is investigated by applying it to solve various test problems. These include: only dispersion effect, DO sag curve, linear variation of cross-sectional area, and application to real flow parameters of Tha Chin river in Thailand. After simplifying the model and considering at steady state flow, certain analytical solutions are possible to derive. We have performed simulations by setting suitable initial and boundary conditions. It can be shown that unsteady solutions of BOD and DO concentrations approach analytical solutions at steady state. This shows the accuracy of our numerical scheme. Furthermore, we can apply the numerical method to study the dynamics of BOD and DO concentrations in the cases of unavailable analytical solution which is typical real flow problems when dealing with arbitrary cross-sectional area or non-constant rate of dispersion that depends the position of channel. Acknowledgments We would like to thank Prof. A.L. Pardhanani at the Earlham College for useful comments and proof-reading throughout this work. References Cockburn, B., Lin, S., Shu, C., 1989. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one-dimensional systems. Journal of Computational Physics 84, 90–113. Jha, R., Bhatia, K.K.S., 2007. On development of refined BOD and DO models for highly polluted Kali river in India. Journal of Environmental Engineering 113, 839–852. Knight, R.L., Payne Jr., V.W.E., Borer, R.E., Clarke, J.R.A., Pries, J.H., 2000. Constructed wetlands for livestock wastewater treatment. Ecological Engineering 15, 41–55. Le, T.V., Roider, E.M., Adrian, D.D., 2007. Laplace transform application to a nontraditional dissolved oxygen model. International Journal of Modeling and Simulation 27, 355–360. Onyejekwe, O.O., Toolsi, S., 2001. Certain aspects of green element computational modelfor BOD–DO interaction. Advances in Water Resources 24, 125–131. Pimpunchat, B., Sweatman, W.L., Wake, G.C., Triampo, W., Parshotam, A., 2009. A mathematical model for pollution in river and its remediation by aeration. Applied Mathematics Letters 22, 304–308. Pollution Control Department (PCD), 2000. Water Quality Standards & Criteria in Thailand, 4th ed. Ministry of Science, Technology and Environment, Thailand. Streeter, H.W., Phelps, E.B., 1925. A Study of the Pollution and Natural Purification of the Ohio river III. Factors Concerned in the Phenomena of Oxidation and Reaeration. U.S. Public Health Service, p. 146. Tyagi, B., Gakkhar, S., Bhargava, D.S., 2009. Modified scheme for one-dimensional BOD–DO models. Journal of Environmental Science and Engineering 51, 145–150. Wang, L.K., Vielkind, D., Wang, M., 1978. Mathematical models of dissolved oxygen concentration in fresh water. Ecological Modelling 5, 115–123. Wu, Z., Li, Z., Chen, G.Q., 2011a. Multi-scale analysis for environmental dispersion in wetland flow. Communications in Nonlinear Science and Numerical Simulation 16, 3168–3178. Wu, Z., Chen, G.Q., Zeng, L., 2011b. Environmental dispersion in a two-zone wetland. Ecological Modelling 222, 456–474.
22
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Montri Maleewong ∗ Sion Hasadsri Department of Mathematics, Faculty of Science, Kasetsart University, Bangkok 10900, Thailand ∗ Corresponding
author at: Department of Mathematics, Faculty of Science, Kasetsart University, Bangkok 10900, Thailand. Tel.: +66 25625555; fax: +66 29428038. E-mail addresses: [email protected], [email protected] (M. Maleewong) 3 April 2012 Available online 9 January 2013
Contents lists available at SciVerse ScienceDirect
Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel
Letter to the Editor a r t i c l e
i n f o
Keywords: BOD and DO Finite element method Dispersion Advection
a b s t r a c t A mathematical model for DO (dissolved oxygen) and BOD (biochemical oxygen demand) in a nonuniform open channel is presented in this work. This model can be used to study transport processes in one-dimensional flow. An accurate numerical scheme based on the finite element method with linear polynomial basis is developed for studying flows both upstream and downstream with respect to a point source. The accuracy of numerical solutions have been studied and compared with certain analytical solutions at steady state. Various test cases including only dispersion effect, DO sag curve, and the linear variation of cross-sectional area are investigated. We have also derived an analytical solution for the last case when the model is simplified under certain conditions. This model has been applied to real flow parameters of Tha Chin river in Thailand. Furthermore, the impacts of flow velocity, kinetic reaction order, and dispersion coefficients have been simulated to observe the dissolved oxygen and biochemical oxygen demand dynamics along the channel. Crown Copyright © 2012 Published by Elsevier B.V. All rights reserved.
1. Introduction Water pollution from human activities including social and economic development, is a crucial problem in environmental management. It directly affects human health and living quality. Also, pollutants from agricultural operations can play a dominant role in the degradation of surface and groundwater quality (Knight et al., 2000). The development of mathematical model enables us to gain better understanding of how to control and predict water quality in rivers or channels. Models generally require various water quality factors to be observed and used as inputs. For example, to model the pollutant distribution in a stream, the important indicators are the level of DO (dissolved oxygen) and BOD (biochemical oxygen demand) concentrations. Flow velocity and saturated oxygen concentration are also the key parameters of the model. The purpose of the model would be to predict the dissolved oxygen concentration and the biochemical oxygen demand concentration along the channel as well as the interaction of these two quantities. Recently, Pimpunchat et al. (2009) have presented a mathematical model of pollution in rivers. They have shown some analytical solutions for various situations of steady state flow. A model involving dissolved oxygen balance in rivers is proposed by Jha and Bhatia (2007), and this model has been improved in various ways. Historically, the famous mathematical model used to predict water quality in rivers was proposed by Streeter and Phelps (1925). They developed the model equations to study how much the declination of dissolved oxygen concentration in a river at certain given points downstream. At steady state flow, the analytical solution can be derived and represented in the form of DO sag curve solution. However, for some extended and complicated models, analytical solutions cannot be obtained. For such cases, efficient numerical methods are required to find approximate solutions of the model. Previously, various numerical methods have been presented for solving pollution models. A finite element method for solving onedimensional problems involving convection-dispersion process and sources of pollutant in river is presented by Onyejekwe and Toolsi (2001). Cockburn et al. (1989) proposed a discontinuous Galerkin method for solving the model equation in conservation law form. They constructed and analyzed a numerical method based on the class of TVB (total variation bounded) schemes. Tyagi et al. (2009) have presented an alternative numerical scheme for predicting the BOD variation in time at various successive distances from a wastewater outfall in a stream. The differential equations include advection, dispersion, and biochemical decay effects. The significance of stream self-purification is proposed by Wang et al. (1978). Water quality is indicated by the level of dissolved oxygen concentration, and the solubility of dissolved oxygen in water. Special emphasis is placed on the modeling of biological reactions and dissolved oxygen saturation concentration in fresh water. The BOD and DO sag models that describe water quality in rivers have been presented by Le et al. (2007). They provide an alternative approach for solving the first-order BOD equation in conjunction with water quality models. In their approach, the differential equation of dissolved oxygen concentration is solved analytically by the application of Laplace transform method. Very recently, multi-scale analysis has been proposed by Wu et al. (2011a) for deriving the analytical solution of multi-scale environmental dispersion in wetland flow for the determination of the dominant variation of contaminant concentration. The environmental dispersion of two-zone shallow wetland extended from single zone is analytically explored by Wu et al. (2011b) in terms of the longitudinal evolution of lateral mean concentration under environmental dispersion. In this work, we present a modified mathematical model for dissolved oxygen concentration interacting with biochemical oxygen demand in a non-uniform open channel flow. The model is in the same form as the equations presented in Pimpunchat et al. (2009). However, in this study, we have considered and developed an accurate numerical method for the case of arbitrary cross-sectional areas on both upstream and downstream flows. The numerical scheme presented is based on the finite element method with linear Lagrange basis functions. The accuracy of numerical solutions depends on the choice of basis function. The advantage of finite element method is that the 0304-3800/$ – see front matter. Crown Copyright © 2012 Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ecolmodel.2012.09.026
12
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Nomenclature P(x,t) X(x,t) Dp Dx l
v A q k1 k2 k ˛ S t n
concentration of BOD (kg m–3 ) concentration of dissolved oxygen (kg m–3 ) dispersion coefficient of pollutant (m2 day–1 ) dispersion coefficient of dissolved oxygen (m2 day–1 ) length of channel (m) average of water velocity (m day–1 ) cross-sectional area (m2 ) rate of pollutant addition (kg m–1 day–1 ) degradation rate coefficient for pollutant (day–1 ) de-aeration rate coefficient for dissolved oxygen (day–1 ) half-saturated oxygen demand concentration for pollutant decay (kg m–3 ) mass transfer of oxygen from air to water (m2 day–1 ) saturated oxygen concentration (kg m–3 ) time (days) number of elements in the finite element method
scheme can be extended directly to higher order by changing the order of basis approximation. In contrast, in the finite difference method derivatives are approximated in terms of difference equations locally and usually suffer from poor approximation of nonlinear term. Using the finite element method, the concentration of dissolved oxygen interacting with biochemical oxygen demand can be approximated at specified positions and time. We have performed various test cases to validate the accuracy of our numerical scheme. Test cases include only dispersion effect, DO sag curve, linear variation of cross-sectional area, and the application to Tha Chin river in Thailand. We have also given an analytical solution for a simplified model. Furthermore, the effects of dispersion, advection, and reaction kinetics are studied by including these terms in the model, and then solving numerically. The mathematical model of dissolved oxygen concentration interacting with biochemical oxygen demand concentration in an open channel is presented in Section 2. The finite element method and time marching scheme are presented in Section 3. Test cases and numerical results are shown in Sections 4 and 5 respectively. Finally, we have made some conclusions in Section 6.
2. Mathematical model We consider a one-dimensional flow along a channel. The spatial variable x is used to describe the distance along the channel. Environmental quantities, such as BOD and DO concentrations, are allowed to vary along the length of the channel. The mathematical model for the interaction between the BOD and DO concentrations along the channel includes the effects of advection, dispersion, and reaction. The model equations can be derived from the conservation of mass for an arbitrary cross-sectional area of the channel. The governing equations can be expressed as follows. ∂2 (AP) X ∂ (vAP) ∂ (AP) = Dp − K1 AP + qH(x), − X +k ∂t ∂x ∂x2
(1)
∂2 (AX) X ∂ (vAX) ∂ (AX) = Dx − K2 AP + ˛(S − X). − X +k ∂t ∂x ∂x2
(2)
This model includes dispersion and advection effects for unsteady processes in arbitrary cross-sectional area which is a function of distance. A reaction term is included in the form of Michaelis–Menten kinetics. The model allows adding pollution which is represented by the Heaviside function H. The reaeration process of dissolved oxygen until reaching the saturated concentration S is represented by the last term in Eq. (2). In this work, we have generalized the model from Pimpunchat et al. (2009) by allowing the cross-sectional area A to be a function of x. In some specific case of flows, analytical solutions are possible to derive. But, for more complex flows, analytical solution cannot be obtained. Thus, we apply the finite element method to solve numerically the model equations. The accuracy of numerical solutions is tested by comparing with some existing analytical solutions for some flow cases.
3. Finite element method The main steps of the finite element formulation for solving Eqs. (1) and (2) are presented as follows.
3.1. Spatial discretization We begin by dividing the computational domain (–l ≤ x ≤l) in n elements, element length is L = 2l/n. The points at x = –l and x = l refer to the end points of upstream and downstream respectively.
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
13
3.2. Approximate function The BOD concentration for each element is approximated by P (x) =
2
i Pi = 1 P1 + 2 P2 = [1 2 ]
i=1
P1
= [] [P] .
P2
(3)
Similarly, the dissolved oxygen concentration is approximated by
X (x) = [1 2 ]
X1
= [] [X] ,
X2
(4)
where 1 = 1 − x/L and 2 = x/L are the linear Lagrange basis functions for each element, [P], and [X] are the values of the unknown functions at each node of the element. 3.3. Weighted residual method The system (1) and (2) is approximated by the standard Galerkin finite element method based on weighted residuals. The system can be written as
L W
∂ (AP) ∂2 (AP) X ∂ (vAP) − Dp + K1 AP − qH(x) + X +k ∂t ∂x ∂x2
dx = 0,
(5)
0
L W
∂ (vAX) ∂ (AX) ∂2 (AX) X − Dx + K2 AP − ˛(S − X) + X +k ∂t ∂x ∂x2
dx = 0,
(6)
0
where W is a weighting function. We choose the weighting function to be the same as the basis function. Thus, the finite element formulation for Eqs. (5) and (6) becomes
L
∂ (AP) i dx − ∂t
L
∂2 (AP) i Dp dx + ∂x2
L
∂ (vAP) dx + i ∂x
L
0
0
0
0
L
L
∂ (AX) i dx − ∂t
∂2 (AX) i Dx dx + ∂x2
0
0
L
X APdx − i K1 X +k
0
i qH(x)dx = 0,
(7)
i ˛(S − X)dx = 0,
(8)
0
L
∂ (vAX) i dx + ∂x
L
X i K2 APdx − X +k
0
L 0
where i refers to the node number in an element (i = 1,2). 3.4. Finite element formulation For the special case when the half-saturated oxygen demand for pollutant decay, k, is negligible (k≈0), Eqs. (7) and (8) simplify to
⎡L ⎤ ⎡ L ⎤ ⎡ L ⎤ dj di dj ⎣ A(x)i j dx⎦ P˙ j + ⎣Dp A(x) dx⎦ Pj + ⎣v A(x)i dx⎦ Pj ⎡
0
+ ⎣K1
dx dx
L
0
⎤
A(x)i j dx⎦ Pj −
0
dx
0
L
,
(9)
,
(10)
qH(x)i dx = 0 0
⎡L ⎤ ⎡ L ⎤ ⎡ L ⎤ dj di dj ⎣ A(x)i j dx⎦ X˙ j + ⎣Dx A(x) dx⎦ Xj + ⎣v A(x)i dx⎦ Xj 0
⎡ + ⎣K2
L
⎤0
dx dx
A(x)i j dx⎦ Pj − ˛S
0
where P˙ j ≈ (Pjn+1 − Pjn )/t, Pj = X˙ j ≈ (Xjn+1 − Xjn )/t, Xj =
L
⎡
L
i dx + ⎣˛
0
Pjn+1 + (1 − )Pjn
0
⎤
dx
i j dx⎦ Xj = 0
0
and
(11)
Xjn+1 + (1 − )Xjn .
(12)
14
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Here Pj and Xj denote the variable’s values at time t = tn and at node xj , t is the time step, and is usually specified in the range of 0 ≤ ≤ 1 and is used to weight the approximate solutions in time. The most commonly used value of is 0 (gives the explicit forward Euler scheme), 0.5 (gives a second-order central implicit method, or Crank-Nicolson-Galerkin), or 1 (leads to the fully implicit method). In an element, we allow the cross-sectional area A to vary as a function of x representing arbitrary cross-sectional area. Thus we approximate A by A = [][A] where [A]T = [A1 A2 ] represents the cross-sectional area at two nodes of the element. The resulting Eqs. (9) and (10) can then be written in matrix forms as
[T ] + t
Ap + [B] + K1 [T ]
Pjn+1 = [T ] + t − 1
Ap + [B] + K1 [T ]
Pjn + t [G] ,
[T ] + t ([Ax ] + [B] + [F2]) Xjn+1 = [T ] + t − 1 ([Ax ] + [B] + [F2]) Xjn − t [[E] − [F1]] ,
where
L
L A(x)i j dx = 12
[T ] =
3A1 + A2
A1 + A2
A1 + A2
A1 + 3A2
0
(13) (14)
,
(15)
and dispersion matrices shown by
L [AP ] == Dp
Dp d dj dx = A(x) i 2L dx dx
0
and,
Dx [Ax ] = 2L
− (A1 + A2 )
− (A1 + A2 )
A1 + A2
A1 + A2
− (A1 + A2 )
− (A1 + A2 )
A1 + A2
The matrix of BOD advection is
L [B] = v
A1 + A2
A(x)i
dj dx
dx =
,
(16)
.
(17)
v
− (2A1 + A2 )
2A1 + A2
6
− (A1 + 2A2 )
A1 + 2A2
0
(18)
The matrix of adding pollution is given by [G] = qH(x)
L 2
1 1
.
(19)
Similarly, matrices [F1], [F2], and [E] are
L i dx = ˛S
[F1] = ˛S 0
1
i j dx =
˛L 3
⎢ ⎣
0
L A(x)i dx = K2
[E] = K2
,
1
⎡
L [F2] = ˛
L 2
0
(20)
1
1 2
1 2
1
L
⎤ ⎥ ⎦,
(21)
xA2 xA1 + A1 − L L
⎡ x⎤
1− L⎦ ⎣ dx.
(22)
x L
0
Finally, after assembling all elements together, we obtain the full form of the matrix system,
[T ]sys + t
sys
+ [B]sys + K1 [T ]sys
sys
+ [B]sys + K1 [T ]sys
Ap
= [T ]sys + t − 1
Ap
= [T ]sys + t − 1
Pjn+1
[T ]sys + t [Ax ]sys + [B]sys + [F2]sys
[Ax ]sys + [B]sys + [F2]sys
(23)
Pjn + t[G]sys ,
Xjn+1
Xjn + t[F1]sys − t[E]sys Pjn+1 − t(1 − )[E]sys Pjn .
(24)
Subscript sys means the ensemble matrices. These two systems can be solved iteratively by the Gauss-Seidel method. For all our computations, the tolerance is set to 10–9 in order to ensure completely converged solutions. The system (23) is solved iteratively before setting the values into the RHS of system (24) to obtain the approximation of Xjn+1 and Pjn+1 at the same time level. When we consider the interaction between BOD and DO, k cannot be neglected. The finite element systems (23) and (24) can be applied, but must be modified slightly. In this case, the Trapezoidal rule is employed for approximating the integral values in matrices [T]sys and [E]sys at time tn .
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
15
3.5. Case of variable dispersion For arbitrary cross-sectional-area as a function of distance, the dispersion coefficient can also be regarded as a function of distance. In this case, the governing equations become, ∂ (AP) ∂ = ∂t ∂x ∂ ∂ (AX) = ∂t ∂x
∂P Dp A ∂x
Dx A
∂X ∂x
−
−
X ∂ (vAP) − K1 AP + qH(x), X +k ∂x
∂ (vAX) X − K2 AP + ˛(S − X). X +k ∂x
Here Dp and Dx are functions of x. To obtain the finite element system, we approximate the dispersion coefficients by
Dp = [1 2 ]
Dp1
and
Dp2
Dx = [1 2 ]
Dx1
,
Dx2
where Dpi and Dxi for i = 1,2 are the dispersion coefficients of pollutant and DO at nodes of the element. The dispersion matrices [Ap ] and [Ax ] in Eqs. (13) and (14) are reformulated as
L [Ap ] =
Dp (x)A(x)
di dj dx, dx dx
Dx (x)A(x)
di dj dx. dx dx
0
L [Ax ] = 0
After integrating, we obtain the modified dispersion matrices for the case of variable dispersion as
Ap
(2Dp1 A1 + Dp2 A2 + Dp2 A1 + 2Dp2 A2 ) = 6L
[Ax ] =
(2Dx1 A1 + Dx2 A2 + Dx2 A1 + 2Dx2 A2 ) 6L
1
−1
−1
1
1
−1
−1
1
,
.
It can be seen that, if Dp1 = Dp2 and Dx1 = Dx2 , it means that Dx and Dp are constants. These resulting matrices are the same as those shown previously in Eqs. (16) and (17). So, for the case of variable dispersion coefficients and variable cross-sectional area, these modified matrices can be applied in place of Eqs. (16) and (17) for flow simulations. 3.6. Convergence of iterative scheme We solve the linear systems (23) and (24) by iterative scheme. So the convergence of our numerical solution depends on the properties of matrix coefficient. If we define, [PO] = [T ]sys + t
Ap
sys
+ [B]sys + K1 [T ]sys ,
and
[DO] = [T ]sys + t [Ax ]sys + [B]sys + [F2]sys ,
the convergence of the linear system (23) and (24) by the Guass–Seidel method depends on the invertibility of matrices [PO] and [DO]. The spectral radius, , which is defined as the largest magnitude of the eigenvalues of the matrix, can be used to analyze the convergence of the iterative scheme. To obtain a converged solution, we require the spectral radius of matrix [PO] denoted by PO , and matrix [DO] denoted by DO less than unity. Deriving explicit forms of stability criteria from these two matrices is difficult and needs to be analyzed in detail further. However, it can be seen that, if the element size is fixed, convergence depends directly on the values of t and , which are the time step and the implicitness of iterative scheme employed in the approximations. 4. Test cases In this section, we present numerical examples for investigating the efficiency of the finite element method. Three test cases are presented: only dispersion effect, the BOD–DO sag-curve behavior, and the linearly varied cross-sectional area. New analytical solutions are derived in this work for checking the accuracy of the numerical method. 4.1. Dispersion effect We consider the case when the dominant effect comes from dispersion that results from different levels of concentration in the longitudinal direction of flow. For simplicity, dispersion coefficient is taken as a constant, and so the exact solution is easily derived. Note that, for more realistic problems, dispersion coefficient should be a function of position, but exact solutions are not available for most problems. / 0 and However, the numerical scheme presented in section 3.5 can be applied for such cases. When k is assumed to be negligible (k≈0), q =
16
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
A is constant, the analytical solutions of BOD and DO concentrations at steady-state on the upstream and downstream flows are presented in Pimpunchat et al. (2009). The solutions are expressed as follows.
⎧ q ı+ˇ ⎪ ı+ˇ)x ( ⎪ 1 − e ,x ≥ 0 ⎪ ⎨ K1 A 2ˇ P(x) = , ⎪ ı+ˇ q ⎪ ı+ˇ)x ( ⎪ e ,x<0 ⎩ K1 A
(25)
2ˇ
v2 +4Dp K1
where ı = 2Dv , ˇ =
2Dp
p
, and
⎧ K2 q ı+ˇ + K2 q ı−ˇ ı + ˇ (ı+ˇ)x ⎪ (+ )x ⎪ + − S − + e − xe ,x ≥ 0 ⎪ ⎨ K1 ˛ K1 2 ˛ 4ˇ A∗ 4ˇ B∗ 2ˇA∗ X(x) = , ⎪ K2 q ı−ˇ ı − ˇ (ı+ˇ)x ı+ˇ − ⎪ (+ )x ⎪ + + e − xe , x < 0 − S ⎩ K1
4ˇ A∗
2 ˛
v2 + 4˛Dx
4ˇ B∗
(26)
2ˇB∗
A , A∗ = 2ADx ı − ˇ − vA , and B∗ = 2ADx ı + ˇ − vA. where = 2Dv x , = 2Dx To compare with the analytical solutions at steady state, we set the initial and boundary conditions in our simulations as
q P(0, t) = K1 A
1−
ı+ˇ 2ˇ
P(x, 0) = 0, X(x, 0) = S,
,
X(0, t) = S −
K2 q K2 q + K1 ˛ K1
+ ı−ˇ ı+ˇ + − 2 ˛ 4ˇ A∗ 4ˇ B∗
,
∂X(l, t) ∂P(−l, t) ∂X(−l, t) ∂P(l, t) = 0, = 0, = 0, and = 0 for t > 0 ∂x ∂x ∂x ∂x
4.2. BOD–DO sag curve When the half-saturated oxygen demand concentration for pollutant decay is negligible (k≈0), and there is no pollution added along the channel, except only the pollutant discharge point at x = 0, the analytical solutions for BOD and DO concentrations on the upstream and downstream flows can be obtained. The solutions at steady state are expressed as follows.
P(x) =
P0 e 2 x , x ≥ 0
(27)
P0 e 1 x , x < 0,
⎧⎛ ⎞ ⎪ ⎪ K K2 P0 P ⎪ 2 0
2 x ⎪ ⎝X0 − S − ⎠ e 4 x + + S, x ≥ 0 ⎪ ˛ ˛e ⎪ 2 2 ⎨ Dx 2 − v 2 − Dx 2 − v 2 − A A ⎛ ⎞ X(x) = ⎪ ⎪ P K K2 P0 ⎪ 2 0
1 x ⎪ ⎝X0 − S − ⎠ e 3 x + + S, x < 0, ⎪ ˛ ˛e ⎪ 2 2 ⎩ Dx 1 − v 1 − Dx 1 − v 1 − A
where
1 =
4 =
v+
v2
+ 4K1 Dp
2Dp
v−
v2 + 4 ˛A Dx 2Dx
,
2 =
v−
A
" v2
+ 4K1 Dp
2Dp
,
3 =
v+
˛ A
v2 + 4 Dx
, and
2Dx
.
To compare numerical results with the analytical solutions at steady state, we set initial and boundary conditions as P(x, 0) = 0, X(x, 0) = S, P0 = P(0, t) = 5, X0 = X(0, t) = S, ∂X(l, t) ∂P(−l, t) ∂X(−l, t) ∂P(l, t) = 0, = 0, = 0, and = 0 for t > 0. ∂x ∂x ∂x ∂x
(28)
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
17
4.3. Linear variation of cross-sectional area In this case, we assume that pollutant is being added along the non-uniform channel at constant rate of q = 0.2 kg m–1 day–1 . The mass transfer of oxygen from air to water ˛ is negligible. The reaction constant is zero (k = 0). To show the ability of our numerical scheme for solving the model in the case of non-uniform cross-sectional area, the function of A is simply assumed to be
# # ⎧ #b − a# ⎪ ⎪ ⎨ a − 200 x, x ∈ [0, 200] A(x) = # # ⎪ #b − a# ⎪ ⎩ a+
200
x, x ∈ [−200, 0) ,
where a and b are the values of cross-sectional area at x = 0 and x = 200 respectively. After solving the governing Eqs. (1) and (2), the analytical solutions can be obtained as follows.
⎧ q q A0 P0 − e 2 x + ⎪ ⎪ K K ⎪ 1 1 ⎪ ,x≥0 ⎨ A P (x) =
⎪ q q ⎪ ⎪ A0 P0 − e 1 x + ⎪ ⎩ K1 K1
X (x) =
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨
A
, x < 0,
q K1 Dx 22 − v 2
K2 A0 P0 − A0 X0 −
(29)
q K1 Dx 22 − v 2
K2 A0 P0 − +
e 2 x
A
,x ≥ 0
⎛ q ⎞ q ⎪ ⎪ K2 A0 P0 − K2 A0 P0 − ⎪ ⎪ K1 ⎠ 5 x K1 ⎪ e + e 1 x ⎪ ⎝A0 X0 − ⎪ 2 − v
2 − v
⎪ D
D
x x 2 2 ⎪ 2 2 ⎪ ⎩ , x < 0,
(30)
A where 5 = Dvx , A0 P0 = A(0,t)P(0,t), and A0 X0 = A(0,t)X(0,t) for t > 0. Initial and boundary conditions in our simulations are given by P(x, 0) = 1, X(x, 0) = S, P0 = P(0, t) = 1, X0 = X(0, t) = S, ∂X(l, t) ∂P(−l, t) ∂X(−l, t) ∂P(l, t) = 0, = 0, = 0, and = 0 for t > 0. ∂x ∂x ∂x ∂x 5. Numerical results In all of our simulations, we set the degradation rate coefficient as K1 = 0.5 day–1 , the de-aeration rate coefficient as K2 = 0.05 day–1 , the saturated oxygen concentration as S = 2 kg m–3 , the water velocity as v = 2 m day–1 , the dispersion coefficient of pollution or BOD as Dp = 5 m2 day–1 , the dispersion coefficient of dissolved oxygen as Dx = 5 m2 day–1 , and the final time of simulations as t = 200 days. This duration is very long, and is used to investigate the convergence of steady solution. By fixing these model parameters, we then simulate the test cases by varying pollutant added along the channel (q), and the mass transfer of oxygen from air to water (␣). 5.1. Dispersion effect We set model parameters to v = 2 m day–1 , q = 0.2 kg m–1 day–1 (for x≥0), q = 0 kg m–1 day–1 (for x<0), Dp = 5 m2 day–1 , Dx = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 . We consider the uniform cross-sectional area in this case. The results of BOD concentration on the upstream and downstream flows at various time steps are shown in Fig. 1. The initial BOD concentration is zero along the channel. The numerical solution of BOD concentration increases dramatically downstream and approaches the analytical solution described in Eq. (25) when k≈50 days. There is no pollution added to the upstream flow, but there is a small dispersion area near the point of x = 0. This shows dispersion effect from the pollution added on downstream. The long-time numerical solutions are in good agreement with the analytical solution at steady state. This shows the accuracy of the present scheme. Moreover, we can apply this numerical method to predict the dynamics of BOD concentration before approaching the steady state in which analytical solutions are not available. The numerical results for the dissolved oxygen concentration are shown in Fig. 2. The initial concentration is given by X(x,0) = S = 2. As time increases, dissolved oxygen concentration decreases exponentially on the downstream channel. Its value, at t≈70days, approaches the analytical solution described in Eq. (26) at steady state. The dissolved oxygen concentration decreases due to the effect of constant pollution adding along the downstream channel. It has no effects to the dissolved oxygen concentration on the upstream flow, except at small area near x = 0. 5.2. BOD–DO sag curve Initially, BOD concentration is zero along the uniform channel, except only at the discharge point of x = 0. At this point, we assume that the BOD concentration is released continuously described by P(0,t) = 5 t > 0. It corresponds to a continuous adding of a point source. There
18
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Fig. 1. BOD concentration versus length when v = 2 m day–1 , q = 0.2 (for x ≥ 0), q = 0 kg m–1 day–1 (for x < 0), Dp = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 .
is no pollutant added on the upstream and downstream flows. The dissolved oxygen concentration on the stream can be increased by the re-aeration process where we set the re-aeration constant as ˛ = 0.01 m2 day–1 . The numerical results of BOD concentration and the analytical solution at steady state are shown in Fig. 3. The unsteady numerical solutions show the dynamics of BOD concentration at various time steps. The BOD concentration gradually increases from the point source and reaches the steady value when t≈9 days. Influent zone downstream (x≈20) is larger than upstream (x≈2) due to the flow direction effect, since we set the flow from left to right. The numerical results of dissolved oxygen concentration and the analytical solution at steady state are shown in Fig. 4. The initial DO concentration is set as X(x,0) = S = 2 which is a constant value along the channel. Numerical simulations show that dissolved oxygen concentration decreases as time increases and approaches the steady value for very large time. DO concentration decreases rapidly in the downstream flow due to the point source effect and then increases again on the far downstream due to re-aeration process. The profiles of DO concentration on the downstream are called the DO sag curve. This simulation shows the effect of pollutant addition and self-recovery area on the downstream. 5.3. Linear variation of cross-sectional area For simplicity in order to derive analytical solution, we assume that the cross-sectional area is varied as a linear function ranging from 10 = m2 to 200 = m2 (for x = 0 to x = l) and from 200 = m2 to 10 = m2 (for x = −l to x = 0). The initial BOD and DO concentrations are set as 1 kg m–3 and 1.6 kg m–3 respectively. Other parameters are given by v = 2 m day–1 , Dp = 5 m2 day–1 , and Dx = 5 m2 day–1 . Pollutant is added on both upstream and downstream by specifying q = 0.2 kg m–1 day–1 , except only at the point of x = 0, where q is zero. There is no re-aeration process in this case, so we set ˛ = 0.
Fig. 2. DO concentration versus length when v = 2 m day–1 , q = 0.2 kg m–1 day–1 (for x ≥ 0), q = 0 kg m–1 day–1 (for x < 0), Dx = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 .
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Fig. 3. BOD concentration versus length when v = 2 m day–1 , Dp = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 .
Fig. 4. DO concentration versus length when v = 2 m day–1 , Dx = 5 m2 day–1 , and ˛ = 0.01 m2 day–1 .
Fig. 5. BOD concentration versus length when v = 2 m day–1 , q = 0.2 kg m–1 day–1 , Dp = 5 m2 day–1 , and ˛ = 0 m2 day–1 .
19
20
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Fig. 6. DO concentration versus length when v = 2 m day–1 , q = 0.2 kg m–1 day–1 , Dx = 2 m2 day–1 , and ˛ = 0 m2 day–1 .
The numerical results of BOD concentration at various time steps are shown in Fig. 5. Due to the constant pollutant loading along the non-uniform channel, the BOD concentration decreases rapidly on the upstream flow near x = 0. The BOD concentrations further decrease and approach zero at the far upstream. Similar results appear downstream. These results reflect the increasing cross-sectional area effects in those positions. The numerical results of dissolved oxygen concentration and the exact solution at steady state from Eq. (30) are shown in Fig. 6. The constant DO is specified continuously at the point of x = 0 meaning the addition of DO concentration. For the same time, constant BOD is also specified at this point. Numerical simulations show that dissolved oxygen concentration decreases due to the constant point source addition both upstream and downstream along the channel. But, in the early downstream flow, the DO concentration does not decrease rapidly like the DO sag curve case. The rate of decrease is smaller due to the cross-sectional area effect that is represented as an increasing function on the left and right of the stream flow. 5.4. Simulation of Tha Chin river In this section, we apply the numerical scheme to simulate the change of DO and BOD concentrations in the Tha Chin river located in Thailand. The variables and flow parameters of this river are reported in Pimpunchat et al. (2009) and PCD (2000). In SI units, the length of river is 325 km, the water velocity is 43,200 m day–1 , the cross-sectional area is 2100 m2 , the rate of pollutant addition along the river is 0.06 kg m–1 day–1 , K1 = 8.27 day–1 , K2 = 44.1 day–1 , ˛ = 16.5 m2 day–1 , and S = 0.01 kg m–3 . We simulate the case of zero dispersion. Graphs of P and DO concentrations along the river are shown in Figs. 7 and 8 respectively. Distance in these figures is shown in units of kilometer, and the unit of concentration is kilogram per hour. We set DO concentration at initial time to 10 mg l–1 . If we add pollutant with constant value along the river, the DO concentration decreases dramatically along the river and approaches the steady solution. For approximately 5 days,
Numerical value of pollution concentration (P) and exact solution at steady state 3500 t = 20 hrs 3000
pollution (P)
2500
2000
1500
1000
500 t = 10 hrs 0
0
50
100
150
200
250
300
350
length of river (x) Fig. 7. BOD concentration versus length of Tha Chin river, dashed line is simulated results and solid line is the exact solution at steady state.
Letter to the Editor / Ecological Modelling 252 (2013) 11–22 6
10
x 10
21
Numerical value of dissolved oxygen (X) and steady solution t = 10 hrs
dissolved oxygen (X)
9.8
9.6
9.4 t = 120 hrs 9.2
9
8.8
0
50
100
150
200
250
300
350
length of river (x) Fig. 8. DO concentration versus length of Tha Chin river, dashed line is simulated results and solid line is the exact solution at steady state.
DO concentration drops from 10 to 9.3 mg l–1 at a distance of 200 km. It can be seen further that DO concentration remains decreased along the downstream river when pollutant added downstream. This shows that our simulator can be used to predict and give some viewpoints of pollution control in rivers. 6. Conclusions In this paper, we present a mathematical model describing the interaction between dissolved oxygen and biochemical oxygen demand in a non-uniform open channel. The model includes both upstream and downstream flows where the separating point can be regarded as a point-source location. Transport processes are assumed to be unsteady in one-dimensional space while flow velocity is constant. This model can also be used to study the dynamics of BOD and DO concentrations when pollution is added along the channel. DO is increased by the re-aeration process until approaching saturated concentration. Thus the solutions from this model provide some viewpoints for ecological management and pollution control in channels. To solve the model, we have developed an accurate scheme based on the finite element method. Numerical solutions in each element are approximated by applying the linear Lagrange basis functions. This then results in a coupled linear system of BOD and DO concentrations at arbitrary partitions of space and time. The accuracy of our numerical scheme is investigated by applying it to solve various test problems. These include: only dispersion effect, DO sag curve, linear variation of cross-sectional area, and application to real flow parameters of Tha Chin river in Thailand. After simplifying the model and considering at steady state flow, certain analytical solutions are possible to derive. We have performed simulations by setting suitable initial and boundary conditions. It can be shown that unsteady solutions of BOD and DO concentrations approach analytical solutions at steady state. This shows the accuracy of our numerical scheme. Furthermore, we can apply the numerical method to study the dynamics of BOD and DO concentrations in the cases of unavailable analytical solution which is typical real flow problems when dealing with arbitrary cross-sectional area or non-constant rate of dispersion that depends the position of channel. Acknowledgments We would like to thank Prof. A.L. Pardhanani at the Earlham College for useful comments and proof-reading throughout this work. References Cockburn, B., Lin, S., Shu, C., 1989. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one-dimensional systems. Journal of Computational Physics 84, 90–113. Jha, R., Bhatia, K.K.S., 2007. On development of refined BOD and DO models for highly polluted Kali river in India. Journal of Environmental Engineering 113, 839–852. Knight, R.L., Payne Jr., V.W.E., Borer, R.E., Clarke, J.R.A., Pries, J.H., 2000. Constructed wetlands for livestock wastewater treatment. Ecological Engineering 15, 41–55. Le, T.V., Roider, E.M., Adrian, D.D., 2007. Laplace transform application to a nontraditional dissolved oxygen model. International Journal of Modeling and Simulation 27, 355–360. Onyejekwe, O.O., Toolsi, S., 2001. Certain aspects of green element computational modelfor BOD–DO interaction. Advances in Water Resources 24, 125–131. Pimpunchat, B., Sweatman, W.L., Wake, G.C., Triampo, W., Parshotam, A., 2009. A mathematical model for pollution in river and its remediation by aeration. Applied Mathematics Letters 22, 304–308. Pollution Control Department (PCD), 2000. Water Quality Standards & Criteria in Thailand, 4th ed. Ministry of Science, Technology and Environment, Thailand. Streeter, H.W., Phelps, E.B., 1925. A Study of the Pollution and Natural Purification of the Ohio river III. Factors Concerned in the Phenomena of Oxidation and Reaeration. U.S. Public Health Service, p. 146. Tyagi, B., Gakkhar, S., Bhargava, D.S., 2009. Modified scheme for one-dimensional BOD–DO models. Journal of Environmental Science and Engineering 51, 145–150. Wang, L.K., Vielkind, D., Wang, M., 1978. Mathematical models of dissolved oxygen concentration in fresh water. Ecological Modelling 5, 115–123. Wu, Z., Li, Z., Chen, G.Q., 2011a. Multi-scale analysis for environmental dispersion in wetland flow. Communications in Nonlinear Science and Numerical Simulation 16, 3168–3178. Wu, Z., Chen, G.Q., Zeng, L., 2011b. Environmental dispersion in a two-zone wetland. Ecological Modelling 222, 456–474.
22
Letter to the Editor / Ecological Modelling 252 (2013) 11–22
Montri Maleewong ∗ Sion Hasadsri Department of Mathematics, Faculty of Science, Kasetsart University, Bangkok 10900, Thailand ∗ Corresponding
author at: Department of Mathematics, Faculty of Science, Kasetsart University, Bangkok 10900, Thailand. Tel.: +66 25625555; fax: +66 29428038. E-mail addresses: [email protected], [email protected] (M. Maleewong) 3 April 2012 Available online 9 January 2013