Skip to main content

A 2-D Barakat-Clark finite difference numerical method and its comparison for water quality simulation modeling

Abstract

Background

This paper provides an advance numerical algorithm to solve both ordinary and partial differential equations of surface water quality models. It uses finite difference methods and structures explicit or implicit or other process forms to solve the water quality model. This study also considers the stability of solutions to obtain more accurate results among those numerical algorithms.

Results

Water quality modeling commonly manifests itself in ordinary and partial differential equations in a realistic world. This study has applied numerical solutions to simulate the changing process of water quality in one and two dimensional spaces or in multiple dimensional spaces. The solutions of these analytical methods are provided in this paper to attest the justifiability of these numerical methods. It demonstrates that the 2-dimensional Barakat-Clark numerical method can be a highly efficient tool in obtaining approximate results of ordinary and partial differential equations, which may prove difficult in finding the accurate solution by using conventional methods. At the same time, the stability analysis corroborated the convergence of those numerical solutions.

Conclusions

This study is the first attempt to compare the multiple numerical methods with the 2-D Barakat-Clark method in the water quality modeling process. The results clearly suggest that the Barakat-Clark method is better in reflecting the accuracy of the water quality modeling with stability for hydrological systems.

Background

In the past few decades, the applications of surface water quality models have been intensely studied as population growth and economic development makes more increases pollution in water resources, particularly in the surface hydrological system. Streeter and Phelps have been progressed the water quality models in 1925. Since then, many scientists (Kellogg 1964; Shampine et al. 1979; Stasa 1985; Christie 1987; Hoffman 1992; Rudi et al. 1997; Cash 2003; Shawgfeh 2004; Walter 2006) have carried out rigorous research in this area. Many water quality models, for example, the BOD model, are formed generally by ordinary or partial differential equations (Davis 1962; Na 1979; Taylor 1982; Evans 1985; Noye Noye 1992; Richard 1997; Moiianty 2004; Liu et al. 2005; Munavalli 2004; Timo et al. 2013). For solving those equations, the traditional methods are using the finite difference method and the finite element method.

In this paper, the application of an explicit and an implicit method to solve the water quality equation includes and tests the first order ordinary differential equation, the second order ordinary differential equation and the second order partial differential equation. To illustrate the reasonability of these solutions, the stability of diffusion equations are also provided in this article. The analytic solutions obtained from these diffusion equations are also compared with other numerical methods. The results suggest that the numerical method can be as good to reflect the accurate solution of water quality diffusion equation in the real world condition.

Methods

Fundamental background of water quality modelling

The transference, diffusion, and degradation mathematical equations of pollutants in the hydrological systems are three-dimensional in form. Therefore, according to the mass balance equation, the initial water quality transfer equation is as follows:

C t + i V i C X i = i X i ϵ T , i C X i ± k S k
(1a)

where

Xi = the three axial X1, X2, X3;

Vi = the velocity of three axial V1, V2, V3;

Sk = the source and sink of contaminants;

C = the concentrate of contaminants;

t = time;

The C, Vi and ϵT,i are the function of Xi. In the left side of the equation (1a), the first item indicates the change rate of concentration of contaminants. The second item on the left side shows the transfer rate of contaminants along the stream axial. The first item, in the right hand side of the equation, presents the dispersion of pollutants along a different axial. There, we are assuming that there is no molecular diffusion and turbulent diffusion of pollutants. The second item in left side indicates the sources of contaminants. The symbol +/− states the contaminants increase or decrease in the hydrological systems. In (1a) is a three-dimensional variable coefficient of the second-order Partial Differential Equation (PDE). The parameters of modeling require much information; especially the variable S on the right side of equation, will cause a serious complication in the solution process. It is difficult to get the analytic a solution. Therefore, people usually simplified this modeling process then used the numerical method to generate the approximate solutions. This paper targets several simplified water quality models and uses finite difference method and structures simulate process forms to solve the problem. It also considers the stability of the solutions to obtain more accurate results among those numerical algorithms.

Numerical methods for water quality modeling

Application of least squares method

If the velocity of the stream is very high, the variable of dispersion will be extremely low, so it can be ignored. Thus, the water shift modeling is presented as follows:

V dC dX = kC
(2a)

The equation (2a) is a first order ordinary differential equation. The analytic solution of it can be obtained as follow:

C = C 0 e k V x
(3a)

Some of the parameters in (3a) can be solved by multiple linear regression method. Concurrently it can use the least squares method to establish the error function. To consider the amount of biochemical oxygen demand (BOD) and disintegration, we use L to indicate BOD, and K 1 present disintegration coefficient of BOD. Then, equation (2a) can is presented as:

dL dt = K 1 L
(4a)

The solution of the above equation is L t = L 0 e K 1 t . Therefore the oxygen consumption is y t = L 0 L t = L 0 1 e k 1 t By using experimental data, applying the least squares method estimates the initial value of L0 and K1. There let K1 = K ′ 1 + h then,

y t = L 0 1 e K 1 + h t = L 0 1 e K 1 t e ht
(5a)

let y(t) = af1 + bf2, and a = L 0 , b = L 0 h , f 1 = 1 e K 1 t , f 2 = t e K 1 t . According to the least squares method to minimize i y t y t testing 2 = i a f 1 + b f 2 y t testing 2 then, we can conclusion those coefficient a and b. The a and b can be indicated as following,

a = f 2 2 f 1 y t f 1 f 2 f 2 y t f 1 2 f 2 2 f 1 f 2 2 , b = f 1 2 f 2 y t f 1 f 2 f 1 y t f 1 2 f 2 2 f 1 f 2 2
(6a)

The L 0 and step-size h can be obtained by the value of a and b. If the h in the equation is larger than 0.001, the result are not stability. Otherwise, the solution is reliable.

Application of explicit and implicit methods

One-dimensional water quality modeling

When the geometry of stream varies greatly, for example the river is very long but it is very shallow; the one-dimensional water quality modeling can be expressed as follows:

E 2 C X 2 V C X + kC = C t
(7a)

Obviously, the function (8a) is a one-dimensional second-order Partial Differential Equation (PDE). It is solved by explicit and implicit methods as following.

Solution of explicit and implicit methods

The equation (8a) can be represented by finite difference. Those equations are expressed as follows:

2 C X 2 = C i + 1 l 2 C i l + C i 1 l Δ X 2
(8a)
C X = C i + 1 l C i 1 l 2 ∆X
(8b)
C t = C i l + 1 C i l ∆t
(8c)

then the explicit method can solve them as follows,

C i l + 1 = E ∆t X 2 + V ∆t 2 ∆X C i 1 l + k∆t 2 E ∆t X 2 + 1 C i l + E ∆t X 2 V ∆t 2 ∆X C i + 1 l
(9a)

the solution of implicit methods can be shown as:

E ∆t X 2 + V ∆t 2 ∆X C i 1 l + 1 + 2 E ∆t X 2 + 1 k∆t C i l + 1 E ∆t X 2 V ∆t 2 ∆X C i + 1 l + 1 = C i l
(10a)

Stability Analysis of explicit and implicit methods

In this section the stability of the equations, (10a) and (11a), are examined. The stability of those methods are tested by using the method of the Fourier stability analysis also known as the von Neumann stability analysis. Let C j l = u l e ijθ , and then the equation (10a) can be written as:

u l + 1 e ijθ = E ∆t X 2 + V ∆t 2 ∆X u l e i j 1 θ + k∆t 2 E ∆t X 2 + 1 u l e ijθ + E ∆t X 2 V ∆t 2 ∆X u l e i j + 1 θ
(11a)

The amplification factor λ equal to:

λ = E ∆t X 2 + V ∆t 2 ∆X ( cos θ i sin θ ) + ( k∆t 2 E ∆t X 2 + 1 ) + ( E ∆t X 2 V ∆t 2 ∆X ) ( cos θ + i sin θ ) = 2 E ∆t X 2 1 cos θ V ∆t ∆X i sin θ + k∆t + 1
(12a)

subject to

λ 2 = k∆t + 1 2 E ∆t X 2 1 cos θ 2 + V ∆t ∆X sin θ 2
(13a)

Let s = 1 − cos θ, s [0, 2] then the equation (13b) can be re-written as:

F s = 4 E 2 t 2 X 2 s 2 2 E ∆t X 2 k∆t + 1 s + k 2 t 2 + 2 k∆t + V 2 t 2 X 2 sin 2 θ
(14a)

Therefore

F s = 8 E 2 t 2 X 2 s 2 E ∆t X 2 k∆t + 1
(15a)
F s = 8 E 2 t 2 X 2 0
(16a)
F s 0 F s 0 F 0 0 F 2 0 k 2 t 2 2 k∆t + V 2 t 2 X 2 sin 2 θ 0 16 E 2 t 2 X 2 4 E ∆t X 2 k∆t + 1 + k 2 t 2 + 2 k∆t + V 2 t 2 X 2 sin 2 θ 0
(17a)

The solution of the explicit is unstable.

On the other hand, the implicit equation (11a) is also tested using the this Fourier analysis as follows. Let C j l = u l e ijθ , the equation can be re-written as:

E ∆t X 2 + V ∆t 2 ∆X u l + 1 e i j 1 θ + 2 E ∆t X 2 + 1 k∆t u l + 1 e ijθ E ∆t X 2 V ∆t 2 ∆X u l + 1 e i j + 1 θ = u l e ijθ

the amplification factor λ is equal to:

λ = 1 / [ E ∆t X 2 + V ∆t 2 ∆X ( cos θ i sin θ ) + ( 2 E ∆t X 2 + 1 k∆t ) ( E ∆t X 2 V ∆t 2 ∆X ) ( cos θ + i sin θ ) ] = 1 / [ 2 E ∆t X 2 1 cos θ + V ∆t ∆X i sin θ k∆t + 1 ]
(18a)

The results show |λ|2 ≤ 1. It indicates that the solution of the implicit method is unconditionally stable. Therefore, through analyzing the stability of equation (12a) and (13a), it is can conclude that the solution of the implicit method will be more accuracate than explicit method in the water quality diffusion equation.

Application the explicit, ADI and Barakat-Clark methods

Two-dimensional water quality modeling

The two-dimensional water quality modeling is represented as follows:

C t + V 1 C X + V 2 C Y = E 1 2 C X 2 + E 2 2 C Y 2 kC
(19a)

As it shows, this model is a second-order Partial Differential Equation (PDE). It will be solved using explicit, ADI and Barakat-Clark methods. The finite difference equation can substituted in the PDE as follows:

C i , j l + 1 C i , j l ∆t / 2 + V 1 C i + 1 , j l C i 1 , j l 2 ∆X + V 2 C i , j + 1 l C i , j 1 l 2 ∆Y = E 1 C i + 1 , j l 2 C i , j l + C i 1 , j l ∆X 2 + E 2 C i , j + 1 l 2 C i , j l + C i , j 1 l ∆Y 2 k C i , j l
(20a)
Solution of explicit, ADI and Barakat-Clark methods

Firstly, equation (20a) is solved by the explicit method. While ∆X is equal to ∆Y, the equation (20a) can be re-written as:

C i , j l = E 1 ∆t 2 X 2 + V 1 ∆t 4 ∆X C i 1 , j l + ( 1 E 1 ∆t 2 X 2 E 2 ∆t 2 X 2 k∆t 2 ) C i , j l + ( E 1 ∆t 2 X 2 V 1 ∆t 4 ∆X ) C i + 1 , j l + E 2 ∆t 2 X 2 + V 2 ∆t 4 ∆X C i , j 1 l + ( E 2 ∆t 2 X 2 V 2 ∆t 4 ∆X ) C i , j + 1 l
(21a)

Secondly, equation (20a) is solved by the ADI method, a two stage time method. At this point, each time step size will be divided into two steps for calculation. The first stage uses a half of step from tl to tl + 1/2 for calculating the result; thus the equation (21a) can be written as:

V 2 ∆t 4 ∆X + E 2 ∆t 2 X 2 C i , j 1 l + 1 / 2 + ( 1 + E 2 ∆t X 2 ) C i , j l + 1 / 2 + ( V 2 ∆t 4 ∆X E 2 ∆t 2 X 2 ) C i , j + 1 l + 1 / 2 = V 1 ∆t 4 ∆X + E 1 ∆t 2 X 2 C i 1 , j l + ( 1 E 1 ∆t X 2 k∆t 2 ) C i , j l + ( E 1 ∆t 2 X 2 V 1 ∆t 4 ∆X ) C i + 1 , j l
(22a)

the second step use anther step size from tl + 1/2 to tl for calculating; the equation (21a) can be written as,

C i , j l + 1 C i , j l + 1 / 2 ∆t / 2 + V 1 C i + 1 , j l + 1 C i 1 , j l + 1 2 ∆X + V 2 C i , j + 1 l + 1 / 2 C i , j 1 l + 1 / 2 2 ∆Y = E 1 C i + 1 , j l + 1 2 C i , j l + 1 + C i 1 , j l + 1 ∆X 2 + E 2 C i , j + 1 l + 1 / 2 2 C i , j l + 1 / 2 + C i , j 1 l + 1 / 2 ∆Y 2 k C i , j l + 1
(23a)

While X is equal to ∆Y, the equation (23a) can be re-written as,

E 1 ∆t 2 X 2 + V 1 ∆t 4 ∆X C i 1 , j l + 1 + ( 1 + E 1 ∆t X 2 + k t 2 ) C i , j l + 1 ( E 1 ∆t 2 X 2 V 1 ∆t 4 ∆X ) C i + 1 , j l + 1 = E 2 ∆t 2 X 2 + V 2 ∆t 4 ∆X C i , j 1 l + 1 / 2 + ( 1 E 2 ∆t X 2 ) C i , j l + 1 / 2 + ( E 2 ∆t 2 X 2 V 2 ∆t 4 ∆X ) C i , j + 1 l + 1 / 2
(24a)

Lastly, equation (20a) is solved using the Barakat-Clark method. According the H. Z. Barakat, the solution of (19a) can be written as (25a) and (25b),

u i , j n + 1 u i , j n ∆t + V 1 u i + 1 , j n u i , j n ∆X + V 2 u i , j + 1 n u i , j n ∆Y = E 1 u i + 1 , j n u i , j n u i , j n + 1 + u i 1 , j n + 1 X 2 + E 2 u i , j + 1 n u i , j n u i , j n + 1 + u i , j 1 n + 1 Δ Y 2 k u i , j n
(25a)
v i , j n + 1 v i , j n ∆t + V 1 v i , j n v i 1 , j n ∆X + V 2 v i , j n v i , j 1 n ∆Y = E 1 v i + 1 , j n + 1 v i , j n + 1 v i , j n + v i 1 , j n X 2 + E 2 v i , j + 1 n + 1 v i , j n + 1 v i , j n + v i , j 1 n Y 2 k v i , j n
(25b)

While ∆X is equal to ∆Y, the equations (25a) and (25b) can be re-written as,

u i , j n + 1 = [ V 1 + V 2 ∆t ∆X ( E 1 + E 2 ) ∆t Δ X 2 k∆t + 1 ) u i , j n + E 1 ∆t X 2 V 1 ∆t ∆X u i + 1 , j n + ( E 2 ∆t X 2 V 2 ∆t ∆X ) u i , j + 1 n + E 1 ∆t X 2 u i 1 , j n + 1 + E 2 ∆t X 2 u i , j 1 n + 1 ] / [ 1 + E + E 2 ∆t X 2 ]
(26a)
v i , j n + 1 = [ V 1 + V 2 ∆t ∆X ( E 1 + E 2 ) ∆t X 2 k t + 1 ) v i , j n + E 1 ∆t X 2 V 1 ∆t ∆X v i 1 , j n + ( E 2 ∆t X 2 V 2 ∆t ∆X ) v i , j 1 n + E 1 ∆t X 2 v i + 1 , j n + 1 + E 2 ∆t X 2 v i , j + 1 n + 1 ] / [ 1 + E 1 + E 2 ∆t X 2 ]
(26b)

The concentration Ci,j of equation (19a) at any time level (n + 1) may be given as,

C i , j n + 1 = u i , j n + 1 + v i , j n + 1 / 2
(27a)

Stability analysis of the explicit, ADI and Barakat-Clark methods

According to the above Fourier stability analysis, the solution of the explicit method (21a) is unstable. For a solution using the ADI method, let C j l = u l e ijθ and substitute it into the equation (24a). This results in:

p = V 2 ∆t 4 ∆X + E 2 ∆t 2 X 2 ( cos θ i sin θ ) + ( 1 + E 2 ∆t X 2 ) + V 2 ∆t 4 ∆X E 2 ∆t 2 X 2 ( cos θ + i sin θ ) = E 2 ∆t ∆X 2 1 cos θ + 1 + V 2 ∆t 2 ∆X i sin θ
(28a)
q = E 1 ∆t 2 X 2 + V 1 ∆t 4 ∆X ( cos θ i sin θ ) + ( 1 E 1 ∆t X 2 k∆t 2 ) + E 1 ∆t 2 X 2 V 1 ∆t 4 ∆X ( cos θ + i sin θ ) = E 1 ∆t ∆X 2 1 cos θ V 1 ∆t 2 ∆X i sin θ + 1 k∆t 2
(28b)

When the amplification factor λ IS equal to λ = q p , s = 1 − cos θ, θ [0, 2] the equations of (28a) and (28b) can be written as:

p = E 2 ∆t X 2 s + 1 + V 2 ∆t 2 ∆X i sin θ
(30a)
q = E 1 ∆t 2 X 2 s V 1 ∆t 2 ∆X i sin θ + 1 k∆t 2
(30b)

By considering |q|2 ≤ |p|2, it is concluded that the solution of the ADI method is stable. The solution stability of the Barakat-Clark method is unconditionally stable because it has been proven by J. A. Clark and H. Z. Barakat (1966).

A case study

The experimental data of 10 groups are showed in the following table for the least squares method. It considers the process which is mentioned in equation (6a) It describes the least squares method to estimate the initial value of BOD and other parameters (Table 1).

Table 1 Experimental data of water samples

Here it assumes that the h = 0.001 and K’1 = 0.1. The results show K1 = 0.3166, L0 = 175.5 by going through the process. The L0 value is the initial value of BOD. Therefore, the water quality model is L(t) = 175.5e− 0.3166t. It is a one-dimensional BOD model.

Now, let’s use those other numerical methods including: explicit, implicit, ADI and Barakat-clark to simulate the degradation of contaminants in a section of stream systems. Under natural conditions, the initial value of concentration of pollutants is assumed to be 0 in this stream section. The inflow from the upstream concentration of contaminants is 150 mg/L. At this point, the equation (8a) is examined. This river section, it is 1000 m long. The water velocity is 17 m/s. The diffusion rate of pollutants E is 360. The degradation rate of the pollutant k is 0.15 in this stream. There is a water quality monitor state in every 200 meters. At the same time, the concentration of pollution is consider to be 0 mg/L at the end of this river, the data is required by the downstream management agencies. Under these conditions, the Department of Environmental Protection wants to know the background concentration of contaminants in this river in order to satisfy the requirement of those downstream cities. Therefore, by using simulation methods, the concentration of pollutants is estimated using the following methods.

The simulation of explicit method

The water quality modeling can be re-written as a standard form of the explicit method.

C i l + 1 = E ∆t X 2 + V ∆t 2 ∆X C i 1 l + k∆t 2 E ∆t X 2 + 1 C i l + E ∆t X 2 V ∆t 2 ∆X C i + 1 l
(31a)

By substituted the initial value, the solution of equations is as follow,

E ∆t X 2 = 0.009 , V ∆t 2 ∆X = 0.0425 , k∆t = 0.15
E ∆t X 2 + V ∆t 2 ∆X = 0.0515 , k t 2 E ∆t X 2 + 1 = 0.832 , E ∆t X 2 V ∆t 2 ∆X = 0.0335

then,

C 1 1 = 0.0515 C 0 0 + 0.832 C 1 0 0.0335 C 2 0 = 7.725
C 2 1 = 0.0515 C 1 0 + 0.832 C 2 0 0.0335 C 3 0 = 0
C 3 1 = 0.0515 C 2 0 + 0.832 C 3 0 0.0335 C 4 0 = 0
C 4 1 = 0.0515 C 3 0 + 0.832 C 4 0 0.0335 C 5 0 = 0
C 1 2 = 0.0515 C 0 1 + 0.832 C 1 1 0.0335 C 2 1 = 14.1522
C 2 2 = 0.0515 C 1 1 + 0.832 C 2 1 0.0335 C 3 1 = 0.3978

The simulation of implicit method

Also, the water quality modeling can be re-written as a standard form of the implicit method.

E ∆t X 2 + V ∆t 2 ∆X C i 1 l + 1 + 2 E ∆t X 2 + 1 + k∆t C i l + 1 E ∆t X 2 V ∆t 2 ∆X C i + 1 l + 1 = C i l
(32a)

Let substitute the value of initial condition. The equations of implicit method is as follow,

E ∆t X 2 + V ∆t 2 ∆X = 0.0515 , 2 E ∆t X 2 + 1 + k∆t = 1.168 , E ∆t X 2 V ∆t 2 ∆X = 0.0335

then

1.168 0.0335 0 0 0.0515 1.168 0.0335 0 0 0.0515 1.168 0.0335 0 0 0.0515 1.168 C 1 1 C 2 1 C 3 1 C 4 1 = 7.725 0 0 0

then these C 1 2 , C 2 2 , C 3 2 , C 4 2 values can be obtained by using above results.

1.168 0.0335 0 0 0.0515 1.168 0.0335 0 0 0.0515 1.168 0.0335 0 0 0.0515 1.168 C 1 2 C 2 2 C 3 2 C 4 2 = 7.725 + C 1 1 C 2 1 C 3 1 C 4 1

The simulation of Barakat-Clark method

The water quality modeling can be re-written as (33a) to (33c), as follows.

u i n + 1 = V ∆t ∆X E ∆t X 2 k∆t + 1 u i n + E ∆t X 2 V ∆t ∆X u i + 1 n + E ∆t X 2 u i 1 n + 1 / 1 + E ∆t Δ X 2
(33a)
1 + E ∆t X 2 v i l + 1 E ∆t X 2 v i + 1 l + 1 = V ∆t ∆X E ∆t X 2 k∆t + 1 v i l E ∆t X 2 + V ∆t ∆X v i 1 l
(33b)
C i , j n + 1 = u i , j n + 1 + v i , j n + 1 / 2
(33c)

The parameters are E ∆t X 2 = 0.009 V ∆t ∆X = 0.0815 k∆t = 0.15 . For u and v; each parameter can be solved as follows:

V ∆t ∆X E ∆t X 2 k∆t + 1 = 0.926 , E ∆t X 2 V ∆t ∆X = 0.076 , 1 + E ∆t X 2 = 1.009 u 1 1 = 0.926 u 1 0 0.076 u 2 0 + 0.009 u 0 1 1.009 = 1.338 u 2 1 = 0.926 u 2 0 0.076 u 3 0 + 0.009 u 1 1 1.009 = 0.012 u 3 1 = 0.926 u 3 0 0.076 u 4 0 + 0.009 u 2 1 1.009 = 0.0001 u 4 1 = 0.926 u 4 0 0.076 u 5 0 + 0.009 u 3 1 1.009 = 0.0000008 u 3 1 = 0.926 u 3 0 0.076 u 4 0 + 0.009 u 2 1 1.009 = 0.0001 u 3 1 = 0.926 u 3 0 0.076 u 4 0 + 0.009 u 2 1 1.009 = 0.0001 u 4 1 = 0.926 u 4 0 0.076 u 5 0 + 0.009 u 3 1 1.009 = 0.0000008
1 + E ∆t X 2 = 1.009 , E ∆t X 2 = 0.009 , V ∆t ∆X E ∆t X 2 k∆t + 1 = 0.756 , E ∆t X 2 V ∆t ∆X = 0.094 1.009 0.009 0 0 0 1.009 0.009 0 0 0 1.009 0.009 0 0 0 1 v 1 1 v 2 1 v 3 1 v 4 1 = 14.1 0 0 0

The simulation of the ADI Method for two-dimensional water quality modeling

The initial condition of the stream has been given as 50 meters long and 40 meters wide. There is no pollution at the beginning. The contaminants came from upstream. This stream section has 4 monitoring points in every 10 meters. The concentrations of pollutants are 150 mg/L at X direction and 150 mg/L at Y direction. The flow velocity is 10 m/s in X direction and 5 m/s in Y direction. The diffusion rate of pollutants E1 and E2 all are equaled to 100. The degradation rate of the pollutant k is 0.15 in this stream. The department of EPA wants to know this section’s water quality every 10 meter long and every 8 meters wide. Therefore, the simulation methods of pollutants in two-dimensional water quality are tested as follows:

According to the ADI method, the two-dimensional water quality equation can be re-written as follows by considering the first half step as (l = 0.5).

V 2 ∆t 2 ∆Y + E 2 ∆t Y 2 C i , j 1 l + 1 / 2 + ( 2 + 2 E 2 ∆t Y 2 ) C i , j l + 1 / 2 + ( V 2 ∆t 2 ∆Y E 2 ∆t Y 2 ) C i , j + 1 l + 1 / 2 = V 1 ∆t 2 ∆X + E 1 ∆t X 2 C i 1 , j l + ( 2 E 1 2 ∆t X 2 k∆t ) C i , j l + ( E 1 ∆t X 2 V 1 ∆t 2 ∆X ) C i + 1 , j l
(34a)

when i = 1, j = 1, 2, 3, 4

1.197 C 1 , 0 0 + 1 / 2 + 3.56 C 1 , 1 0 + 1 / 2 0.363 C 1 , 2 0 + 1 / 2 = 0.75 C 0 , 1 0 + 0.925 C 1 , 1 0 + 0.25 C 2 , 1 0
1.197 C 1 , 1 0 + 1 / 2 + 3.56 C 1 , 2 0 + 1 / 2 0.363 C 1 , 3 0 + 1 / 2 = 0.75 C 0 , 2 0 + 0.925 C 1 , 2 0 + 0.25 C 2 , 2 0
1.197 C 1 , 2 0 + 1 / 2 + 3.56 C 1 , 3 0 + 1 / 2 0.363 C 1 , 4 0 + 1 / 2 = 0.75 C 0 , 3 0 + 0.925 C 1 , 3 0 + 0.25 C 2 , 3 0
1.197 C 1 , 3 0 + 1 / 2 + 3.56 C 1 , 4 0 + 1 / 2 0.363 C 1 , 5 0 + 1 / 2 = 0.75 C 0 , 4 0 + 0.925 C 1 , 4 0 + 0.25 C 2 , 4 0

There C 1 , 1 0 + 1 / 2 . C 1 , 2 0 + 1 / 2 , C 1 , 3 0 + 1 / 2 , C 1 , 4 0 + 1 / 2 , they can be obtained by above equations. Then, it substituted i = 2, j = 1, 2, 3, 4, the equation (34a) can be written as following,

1.197 C 2 , 0 0 + 1 / 2 + 3.56 C 2 , 1 0 + 1 / 2 0.363 C 2 , 2 0 + 1 / 2 = 0.75 C 1 , 1 0 + 0.925 C 2 , 1 0 + 0.25 C 3 , 1 0
1.197 C 2 , 1 0 + 1 / 2 + 3.56 C 2 , 2 0 + 1 / 2 0.363 C 2 , 3 0 + 1 / 2 = 0.75 C 1 , 2 0 + 0.925 C 2 , 2 0 + 0.25 C 3 , 2 0
1.197 C 2 , 2 0 + 1 / 2 + 3.56 C 2 , 3 0 + 1 / 2 0.363 C 2 , 4 0 + 1 / 2 = 0.75 C 1 , 3 0 + 0.925 C 2 , 3 0 + 0.25 C 3 , 3 0
1.197 C 2 , 3 0 + 1 / 2 + 3.56 C 2 , 4 0 + 1 / 2 0.363 C 2 , 5 0 + 1 / 2 = 0.75 C 1 , 4 0 + 0.925 C 2 , 4 0 + 0.25 C 3 , 4 0
C 2 , 1 0 + 1 / 2 . C 2 , 2 0 + 1 / 2 , C 2 , 3 0 + 1 / 2 , C 2 , 4 0 + 1 / 2

, they can be obtained by above equations. Similarly, let substitute i = 3, 4, and j = 1, 2, 3, 4.

The second half of the ADI method is based on the results of the first step. By using i = 1, j = 1, 2, 3, 4, the equation (34a) can be written as,

0.0165 C 0 , 1 1 + 2.218 C 1 , 1 1 0.0015 C 2 , 1 1 = 50.5 C 1 , 0 0 + 1 / 2 98 C 1 , 1 0 + 1 / 2 49.5 C 1 , 2 0 + 1 / 2
0.0165 C 1 , 1 1 + 2.218 C 2 , 1 1 0.0015 C 3 , 1 1 = 50.5 C 2 , 0 0 + 1 / 2 98 C 2 , 1 0 + 1 / 2 49.5 C 2 , 2 0 + 1 / 2
0.0165 C 2 , 1 1 + 2.218 C 3 , 1 1 0.0015 C 4 , 1 1 = 50.5 C 3 , 0 0 + 1 / 2 98 C 3 , 1 0 + 1 / 2 49.5 C 3 , 2 0 + 1 / 2
0.0165 C 3 , 1 1 + 2.218 C 4 , 1 1 0.0015 C 5 , 1 1 = 50.5 C 4 , 0 0 + 1 / 2 98 C 4 , 1 0 + 1 / 2 49.5 C 4 , 2 0 + 1 / 2

These C 1 , 1 1 . C 2 , 1 1 , C 3 , 1 1 , C 4 , 1 1 can be solved using above equations. Substituting j = 1, i = 1, 2, 3, 4, the equation (34a) can be written as,

E 1 ∆t X 2 + V 1 ∆t ∆X C i 1 , j l + 1 + ( 2 + E 1 2 ∆t X 2 + k∆t ) C i , j l + 1 ( E 1 ∆t X 2 V 1 ∆t 2 ∆X ) C i + 1 , j l + 1 = E 2 ∆t Y 2 + V 2 ∆t 2 ∆Y C i , j 1 l + 1 / 2 + ( 2 E 2 ∆t Y 2 ) C i , j l + 1 / 2 + ( E 2 ∆t Y 2 V 2 ∆t 2 ∆Y ) C i , j + 1 l + 1 / 2
(34b)
0.75 C 0 , 1 1 + 3.075 C 1 , 1 1 0.25 C 2 , 1 1 = 1.197 C 1 , 0 0 + 1 / 2 + 1.22 C 1 , 1 0 + 1 / 2 + 0.363 C 1 , 2 0 + 1 / 2
0.75 C 1 , 1 1 + 3.075 C 2 , 1 1 0.25 C 3 , 1 1 = 1.197 C 2 , 0 0 + 1 / 2 + 1.22 C 2 , 1 0 + 1 / 2 + 0.363 C 2 , 2 0 + 1 / 2
0.75 C 2 , 1 1 + 3.075 C 3 , 1 1 0.25 C 4 , 1 1 = 1.197 C 3 , 0 0 + 1 / 2 + 1.22 C 3 , 1 0 + 1 / 2 + 0.363 C 3 , 2 0 + 1 / 2
0.75 C 3 , 1 1 + 3.075 C 4 , 1 1 0.25 C 5 , 1 1 = 1.197 C 4 , 0 0 + 1 / 2 + 1.22 C 4 , 1 0 + 1 / 2 + 0.363 C 4 , 2 0 + 1 / 2

Let substitute j = 2, i = 1, 2, 3, 4, the equation (34a) can be written as,

0.75 C 0 , 2 1 + 3.075 C 1 , 2 1 0.25 C 2 , 2 1 = 1.197 C 1 , 1 0 + 1 / 2 + 1.22 C 1 , 2 0 + 1 / 2 + 0.363 C 1 , 3 0 + 1 / 2
0.75 C 1 , 2 1 + 3.075 C 2 , 2 1 0.25 C 3 , 2 1 = 1.197 C 2 , 1 0 + 1 / 2 + 1.22 C 2 , 2 0 + 1 / 2 + 0.363 C 2 , 3 0 + 1 / 2
0.75 C 2 , 2 1 + 3.075 C 3 , 2 1 0.25 C 4 , 2 1 = 1.197 C 3 , 1 0 + 1 / 2 + 1.22 C 3 , 2 0 + 1 / 2 + 0.363 C 3 , 3 0 + 1 / 2
0.75 C 3 , 2 1 + 3.075 C 4 , 2 1 0.25 C 5 , 2 1 = 1.197 C 4 , 1 0 + 1 / 2 + 1.22 C 4 , 2 0 + 1 / 2 + 0.363 C 4 , 3 0 + 1 / 2

Thus, C 1 , 2 1 . C 2 , 2 1 , C 3 , 2 1 , C 4 , 2 1 can be obtained by above equations. Similarly, the other C value is calculated by substituting i =3, 4.

The simulation of Barakat-Clark method for two-dimensional modeling

This model is similar to the one-dimensional water quality modeling. If it considers the two different directions (x,y). The group of equations can be written as,

u i , j l + 1 = { [ V 1 ∆t ∆X + V 2 ∆t ∆Y E 1 ∆t Δ X 2 E 2 ∆t Y 2 k t + 1 u i , j l + ( E 1 ∆t X 2 V 1 ∆t ∆X ) u i + 1 , j l + E 2 ∆t Y 2 V 2 ∆t ∆Y u i , j + 1 l + E 1 ∆t X 2 u i 1 , j l + 1 + E 2 ∆t Y 2 u i , j 1 l + 1 ] } / [ 1 + E 1 ∆t X 2 + E 2 ∆t Y 2 ]
(35a)

in this case, the u values are as follows

u 1 , 1 1 = 0.312 u 1 , 1 0 + 0.25 u 2 , 1 0 + 0.363 u 1 , 2 0 + 0.5 u 0 , 1 1 + 0.78 u 1 , 0 1 / 2.28 = 84.21
u 1 , 2 1 = 0.312 u 1 , 2 0 + 0.25 u 2 , 2 0 + 0.363 u 1 , 3 0 + 0.5 u 0 , 2 1 + 0.78 u 1 , 1 1 / 2.28 = 61.7
u 1 , 3 1 = 0.312 u 1 , 3 0 + 0.25 u 2 , 3 0 + 0.363 u 1 , 4 0 + 0.5 u 0 , 3 1 + 0.78 u 1 , 2 1 / 2.28 = 54
u 1 , 4 1 = 0.312 u 1 , 4 0 + 0.25 u 2 , 4 0 + 0.363 u 1 , 5 0 + 0.5 u 0 , 4 1 + 0.78 u 1 , 3 1 / 2.28 = 51.37
u 2 , 1 1 = 0.312 u 2 , 1 0 + 0.25 u 3 , 1 0 + 0.363 u 2 , 2 0 + 0.5 u 1 , 1 1 + 0.78 u 2 , 0 1 / 2.28 = 69.783
u 2 , 2 1 = 0.312 u 2 , 2 0 + 0.25 u 3 , 2 0 + 0.363 u 2 , 3 0 + 0.5 u 1 , 2 1 + 0.78 u 2 , 1 1 / 2.28 = 37.4
u 2 , 3 1 = 0.312 u 2 , 3 0 + 0.25 u 3 , 3 0 + 0.363 u 2 , 4 0 + 0.5 u 1 , 3 1 + 0.78 u 2 , 2 1 / 2.28 = 24.64
u 2 , 4 1 = 0.312 u 2 , 4 0 + 0.25 u 3 , 4 0 + 0.363 u 2 , 5 0 + 0.5 u 1 , 4 1 + 0.78 u 2 , 3 1 / 2.28 = 19.695

The rest of u values can be calculated by substituting i =3, 4.

E 1 ∆t X 2 v i + 1 , j l + 1 + E 2 ∆t Y 2 v i , j + 1 l + 1 1 + E 1 ∆t X 2 + E 2 ∆t Y 2 v i , j l + 1 = V 1 ∆t ∆X + V 2 ∆t ∆Y + E 1 ∆t X 2 + E 2 ∆t Y 2 + k∆t 1 v i , j l ( V 1 ∆t ∆X + E 1 ∆t X 2 ) v i 1 , j l V 2 ∆t ∆Y + E 2 ∆t Δ Y 2 v i , j 1 l
(35b)

The v values are calculated by the following process:

0.5 v 2 , 1 1 + 0.78 v 1 , 2 1 1.022 v 1 , 1 1 = 1.022 v 1 , 1 0 0.75 v 0 , 1 0 1.197 v 1 , 0 0 = 292.05
0.5 v 2 , 2 1 + 0.78 v 1 , 3 1 1.022 v 1 , 2 1 = 1.022 v 1 , 2 0 0.75 v 0 , 2 0 1.197 v 1 , 1 0 = 112.5
0.5 v 2 , 3 1 + 0.78 v 1 , 4 1 1.022 v 1 , 3 1 = 1.022 v 1 , 3 0 0.75 v 0 , 3 0 1.197 v 1 , 2 0 = 112.5
0.5 v 2 , 4 1 + 0.78 v 1 , 5 1 1.022 v 1 , 4 1 = 1.022 v 1 , 4 0 0.75 v 0 , 4 0 1.197 v 1 , 3 0 = 112.5
0.5 v 3 , 1 1 + 0.78 v 2 , 2 1 1.022 v 2 , 1 1 = 1.022 v 2 , 1 0 0.75 v 1 , 1 0 1.197 v 2 , 0 0 = 179.55
0.5 v 3 , 2 1 + 0.78 v 2 , 3 1 1.022 v 2 , 2 1 = 1.022 v 2 , 2 0 0.75 v 1 , 2 0 1.197 v 2 , 1 0 = 112.5
0.5 v 3 , 3 1 + 0.78 v 2 , 4 1 1.022 v 2 , 3 1 = 1.022 v 2 , 3 0 0.75 v 1 , 3 0 1.197 v 2 , 2 0 = 0
0.5 v 3 , 4 1 + 0.78 v 2 , 5 1 1.022 v 2 , 4 1 = 1.022 v 2 , 4 0 0.75 v 1 , 4 0 1.197 v 2 , 3 0 = 0

When these values of u and v were calculated, the final solution of the Barakart-Clark method can be obtained by using average value of those two matrixes.

Results and discussion

Figure 1 shows the solutions of explicit, implicit and Barakat-Clark methods. It proves that the numerical solution obtained by the Barakat-Clark method approached the real world process even through the explicit and implicit results are stable in this case. Figure 2 shows the solutions obtained using the ADI and Barakat-Clark methods for the 2-D water quality modelling process. In Figure 2a, it is the results of ADI methods in three different times. Similarly, Figure 2b shows the results of the Barakat-Clark method. The ADI method clearly shows that the results are not as stable as the results of Barakat-Clark method. In fact, during the final time, the ADI result is beyond of boundary conditions. However, it may have been caused by a large space of x or y. However, under same steps and space of x and y, the Barakat-Clark method obtained the better accuracy results. The Barakat-Clark method uses less computing time to get those results. This method has the advantage of solving the three dimensions time-dependent equations which is not convenient for ADI method.

Figure 1
figure 1

The results of explicit, implicit and barakat-clark methods for 1-D water quality modeling.

Figure 2
figure 2

The ADI results of 2-D water quality simulation at begin, middle and final times. The Barakak-Clark results of 2-D water quality simulation at begin, middle and final times.

The data results (Tables 2, 3, 4, 5 and 6) of the above equations show the detailed numerical solutions of the water quality modelling process. The findings in Table 4 show smoother curves of water quality because of different method. Although the explicit method obtained stable results in this case, when X value is changed to a smaller value, the results will be extremely unstable. It means that the stability of the solution is depends on the size of grid for explicit method. It is obvious to see that the solution of the Barakat-Clark method is more stable, even though the size of the grid was changed. This is because the Barakat-Clark numerical method has been considered average values or mean values of implicit and explicit data. Therefore, the Barakat-Clark methodology is more accurate than others in the case of the water quality simulation process.

Table 2 Concentration results of pollutants using the explicit method for 1-D water quality modeling
Table 3 Concentration results of pollutants using the implicit method for 1-D water quality modeling
Table 4 Concentration results of pollutants using the Barakat-Clark method for 1-D water quality modeling
Table 5 Results of ADI method for 2-D water quality modeling in final time ( x/y direction)
Table 6 Results of Barakat-Clark method for 2-D water quality modeling in final time ( x/y direction)

Conclusions

In this study, the finite difference numerical methods are applied to the water quality modeling processes. They simulate the concentration process of pollutants in hydrological systems. As one of the existing numerical methods, the Barakat-Clark has obtained higher accuracy results with stability than the other methods. By comparing the solution of Barakat-Clark method with the other numerical methods, its results better reflect the reality of the simulation without oversimplify solving the process. The water quality modeling process, which estimates the real world conditions, can be effectively solved by means of the Barakat-Clark method. It enables the managers of water authorities to know the concentration of pollution in surface water systems and to conveniently use those solutions as a reference in their hydrological systems. Consequently, the accurate of water quality simulation processes can be enhanced. The results of the case study indicates that the solution of Barakat-Clark method has proved that this method is better by comparison than the others for water quality modeling. Solutions of the numerical method actually reflect a compromise between the finite difference method and requirement of accuracy.

In the real world, administrators of water management agencies are paying a high price for water quality to guarantee the safe conditions of environmental water. This means that authorities desire to accurate monitoring data to determine the stream or river situations. This may cause a higher economic cost. However, simulation values of water quality are reliable within the Barakat-Clark method. Its advantages have been clearly shown in the above case study. It shows that the results of the Barakat-Clark method is more reliable and more accurate with process stability. On the other hand, other numerical methods such as ADI method produced a different results under same circumstances. Therefore, the Barakat-Clark method can be considered as a better finite method in the 2-D water modeling system, and this paper is the first attempt to compare the Barakat-Clark 2-D method with other multiple numerical algorithms in the water quality modeling process. From the above applications, it is obvious that the Barakat-Clark method shares higher stability results with the same environmental conditions. However, in this study, the boundary conditions have been pre-defined. Therefore, further studies of this method have to correspond to the variation of boundary conditions in real world. Finally, this the 2-D Barakat-Clark method well reflects the accuracy of simulation process, and thus, it can provide an exact reference for water quality modeling in hydrological systems.

References

  • Barakat HZ, Clark JA: On the solution of the diffusion equations by numerical methods. J Heat Tran 1966, 11: 421–427.

    Article  Google Scholar 

  • Cash JR: Efficient numerical methods for the solution of stiff initial-value problems and differential algebraic equations. Proc R Soc Lond 2003, 459: 797–815. 10.1098/rspa.2003.1130

    Article  Google Scholar 

  • Christie MA, Bond DJ: Detailed simulation of unstable processes in miscible flooding. Reserv Eng 1987, 2: 514–522.

    Article  Google Scholar 

  • Davis HT: Induction to nonlinear differential and integral equations. New York: Dover; 1962.

    Google Scholar 

  • Evans DJ, Abdullah AR: A new explicit method for diffusion-convection equation. Comp Maths Applis 1985, 11: 145–154. 10.1016/0898-1221(85)90143-9

    Article  Google Scholar 

  • Hoffman J: Numerical methods for engineers and scientists. New York: McGraw-Hill; 1992.

    Google Scholar 

  • Kellogg RB: An alternating direction method for operator equations. SIAM 1964, 12: 848–854.

    Google Scholar 

  • Liu WC, Jantai K, Albert Y: Modelling hydrodynamics and water quality in the separation waterway of the Yulin offshore industrial park, Tainwan. Environ Model Softw 2005, 20: 309–328. 10.1016/j.envsoft.2003.12.013

    Article  Google Scholar 

  • Moiianty RK: An unconditionally stable difference scheme for the one-space-dimensional linear hyperbolic equation. Appl Math Lett Letters 2004, 17: 101–105. 10.1016/S0893-9659(04)90019-5

    Article  Google Scholar 

  • Munavalli GR, Mohan Kumar MS: Modified Lagrangian method for modeling water quality in distribution systems. Water Res 2004, 38: 973–2988.

    Google Scholar 

  • Na TY: Computational methods in engineering boundary value problem. New York: Academic Press; 1979.

    Google Scholar 

  • Noye BJ, Hayman KJ: Explicit two-level finite difference methods for the two-dimensional diffusion equation. J Comput Math 1992, 42: 223–236.

    Google Scholar 

  • Richard JB, Robert S: Expanded stability through higher accuracy for time-centered advection schemes. Springer, Herdelberg: American Meteorological Society; 1997:6.

    Google Scholar 

  • Rudi R, Matjaz C: Hydrodynamic and water quality modeling: An experience. Ecol Model 1997, 101: 195–207. 10.1016/S0304-3800(97)00047-1

    Article  Google Scholar 

  • Shampine LF, Gear CW: A user’s view of solving stiff ordinary differential equations. SIAM Rev 1979, 21: 1. 10.1137/1021001

    Article  Google Scholar 

  • Shawgfeh N, Kaya D: Comparing numerical methods for the solutions of systems of ordinary differential equations. Appl Math Lett 2004, 17: 323–328. 10.1016/S0893-9659(04)90070-5

    Article  Google Scholar 

  • Stasa FL: Applied finite element analysis for engineers. New York: Holt Rinehart and Winston; 1985.

    Google Scholar 

  • Taylor JR: An introduction to error analysis. CA: Mill Valley; 1982.

    Google Scholar 

  • Timo T, Sampsa K, Ville K, Matthieu M, Peng CY: Water quality analysis using an inexpensive device and a mobile phone. Env Syst Res 2013, 2: 9. 10.1186/2193-2697-2-9

    Article  Google Scholar 

  • Walter MG: Use of distribution system water quality models in support of water security. Springer, Heidelerg: Security of Water Supply System; 2006:39–50.

    Google Scholar 

Download references

Acknowledgements

This research was supported by the Program for Innovative Research Team (IRT1127), the Natural Science and Engineering Research Council of Canada, the MOE Key Project Program (311013), Natural Science Foundation of China (51109181) and the Scientific Research fund project (XKJJ201105) of Xiamen University of Technology. The writers are very grateful to the editor and the anonymous reviewers for their insightful comments and suggestions.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Lei Jin.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

Author LJ composed this paper, and others have revised it many times for publication. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Jin, L., Yang, A., Wang, L. et al. A 2-D Barakat-Clark finite difference numerical method and its comparison for water quality simulation modeling. Environ Syst Res 2, 11 (2013). https://doi.org/10.1186/2193-2697-2-11

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2193-2697-2-11

Keywords