# A Numerical Treatment of Nondimensional Form of Water Quality Model in a Nonuniform Flow Stream Using Saulyev Scheme

Authors: Nopparat Pochai
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

The stream water quality model of water quality assessment problems often
involves numerical methods to solve the equations. The governing equation
of the uniform flow model is one-dimensional advection-dispersion-reaction
equations (ADREs). In this paper, a better finite difference scheme for solving
ADRE is focused, and the effect of nonuniform water flows in a stream is
considered. Two mathematical models are used to simulate pollution due to
sewage effluent. The first is a hydrodynamic model that provides the velocity
field and elevation of the water flow. The second is a advection-dispersion-reaction
model that gives the pollutant concentration fields after input of the
velocity data from the hydrodynamic model. For numerical techniques, we
used the Crank-Nicolson method for system of a hydrodynamic model and
the explicit schemes to the dispersion model. The revised explicit schemes are
modified from two computation techniques of uniform flow stream problems:
forward time central space (FTCS) and Saulyev schemes for dispersion model.
A comparison of both schemes regarding stability aspect is provided so as to
illustrate their applicability to the real-world problem.


---

## Body

## 1. Introduction

Field measurement and mathematical simulation are methods to detect the amount of pollutant in water area. For the shallow water mass transport problems that presented in [1], the method of characteristics has been reported as being applied with success, but it presents in real cases some difficulties. In [2], the finite element method for solving a steady water pollution control to achieve a minimum cost is presented. The numerical techniques for solving the uniform flow of stream water quality model, especially the one-dimensional advection-diffusion-reaction equation, are presented in [3–7]. A two-dimensional model for natural convection in shallow water that reduces to a degenerated elliptic equation for the pressure, an explicit formula for horizontal components of the velocity and a vertical diffusion for the vertical component, is derived [8]. In [9], a rigorous nonlinear mathematical model is used to explain the seasonal variability of plankton in previous shallow coastal lagoons. The particle trajectories in a constant vorticity shallow water term flow over a flat bed as periodic waves propagate on the water’s free surface are investigated in [10].The most of nonuniform flow model requires data concerned with velocity of the current at any point and any time in the domain. The hydrodynamics model provides the velocity field and tidal elevation of the water. In [11–13], they used the hydrodynamics model and convection-diffusion equation to approximate the velocity of the water current in a bay, a uniform reservoir, and a channel, respectively.The numerical techniques to solve the nonuniform flow of stream water quality model, one-dimensional advection-diffusion-reaction equation, is presented in [14] by using the fully implicit schemes: Crank-Nicolson method system of hydrodynamic model and backward time central space (BTCS) for dispersion model.The finite difference methods, including both explicit and implicit schemes, are mostly used for one-dimensional problems such as in longitudinal river systems [15]. Researches on finite difference schemes have considered on numerical accuracy and stability. There are several high quality numerical schemes, such as QUICK/QUICKest schemes, Lax-Wendroff scheme, Crandall scheme, and Dufort-Frankel scheme have been developed to enhance model performances. These schemes have outstanding stability and high-order accuracy. They are requirements for advection-dominated systems. Although these schemes need boundary and initial conditions that make them difficult to use. They need more computing effort since iterations for more grids are involved in each computation step. For example, the QUICKest scheme uses a three-point upstream-weighted quadratic interpolation and needs the stop criteria controlled iterations for each grid in order to enhance accuracy. The scheme carries out a heavy computing load. Since it involves two upstream points, the upper boundary conditions need to be defined carefully before starting computation [4].The simple finite difference schemes become more inviting for general model use. The simple explicit schemes include Forward-Time/Centered-Space (FTCS) scheme and the Saulyev scheme. These schemes are either first-order or second-order accurate [4] and have the advantages of simplicity in coding and time effectiveness in computing without losing too much accuracy and thus are preceding for several model applications.In this paper, we will use more economical computation techniques than the method in [14]. For numerical techniques, we used the Crank-Nicolson method to the system of hydrodynamic model and the explicit schemes to the dispersion model. The revised explicit schemes are modified from two computation techniques of uniform flow stream problems: forward in time/central in space (FTCS) and Saulyev schemes.The results from hydrodynamic model are data of the water flow velocity for advection-diffusion-reaction equation which provides the pollutant concentration field. The term of friction forces due to the drag of sides of the stream is considered. The theoretical solution of the model at the end point of the domain that guaranteed the accurate of the approximate solution is presented in [13, 14].The stream has a simple one-space dimension as shown in Figure1. Averaging the equation over the depth, discarding the term due to Coriolis force, it follows that the one-dimensional shallow water and advection-diffusion-reaction equations are applicable. We use the Crank-Nicolson method and the forward in time/central in space (FTCS) and Saulyev schemes to approximate the velocity and the tidal elevation and the concentration from the foresaid models, respectively.Figure 1
The shallow water system.


## 2. Mathematical Model

### 2.1. The Hydrodynamic Model

The continuity and momentum equations are governed by the hydrodynamic behavior of the stream. If we average the equations over the depth, discarding the term due to Coriolis parameter, shearing stresses, and surface wind [11, 13, 14, 16], we obtain the one-dimensional shallow water equations: (2.1)∂ζ∂t+∂∂x[(h+ζ)u]=0,∂u∂t+g∂ζ∂x=0,
where x is longitudinal distance along the stream (m), t is time (s), h(x) is the depth measured from the mean water level to the stream bed (m), ζ(x,t) is the elevation from the mean water level to the temporary water surface or the tidal elevation (m/s), and u(x,t) is the velocity components (m/s), for all x∈[0,l].Assume thath is a constant and ζ≪h. Then (2.1) leads to (2.2)∂ζ∂t+h∂u∂x=̇0,∂u∂t+g∂ζ∂x=0.
We will consider the equation in dimensionless problem by letting U=u/gh, X=x/l, Z=ζ/h, and T=tgh/l. Substituting them into (2.2) leads to (2.3)∂Z∂T+∂U∂X=0,∂U∂T+∂Z∂X=0.
We now introduce a damping term into (2.3) to represent frictional forces due to the drag of sides of the stream, thus (2.4)∂Z∂T+∂U∂X=0,∂U∂T+∂Z∂X=-U,
with the initial conditions at t=0 and 0≤X≤1 being specified: Z=0 and U=0. The boundary conditions for t>0 are specified: Z=eit at X=0 and ∂Z/∂X=0 at X=1. The system of (2.4) is called the damped equation. We solve the damped equation by using the finite difference method. In order to solve (2.4) in [0,1]×[0,T], for convenient using u,d for U and Z, respectively, (2.5)∂u∂t+∂d∂x=-u,∂d∂t+∂u∂x=0,
with the initial conditions u=0, d=0 at t=0 and the boundary conditions d(0,t)=f(t) and ∂d/∂x=0 at x=1.


### 2.2. Dispersion Model

In a stream water quality model, the governing equations are the dynamic one-dimensional advection-dispersion-reaction equations (ADREs). A simplified representation by averaging the equation over the depth is shown in [3–5, 7, 14] as (2.6)∂C∂t+u∂C∂x=D∂2C∂x2-KC,
where C(x,t) is the concentration averaged in depth at the point x and at time t, D is the diffusion coefficient, K is the mass decay rate, and u(x,t) is the velocity component, for all x∈[0,1]. We will consider the model with following conditions. The initial condition C(x,0)=0 at t=0 for all x>0. The boundary conditions C(0,t)=C0 at x=0 and ∂C/∂x=0 at x=1, where C0 is a constant.


## 2.1. The Hydrodynamic Model

The continuity and momentum equations are governed by the hydrodynamic behavior of the stream. If we average the equations over the depth, discarding the term due to Coriolis parameter, shearing stresses, and surface wind [11, 13, 14, 16], we obtain the one-dimensional shallow water equations: (2.1)∂ζ∂t+∂∂x[(h+ζ)u]=0,∂u∂t+g∂ζ∂x=0,
where x is longitudinal distance along the stream (m), t is time (s), h(x) is the depth measured from the mean water level to the stream bed (m), ζ(x,t) is the elevation from the mean water level to the temporary water surface or the tidal elevation (m/s), and u(x,t) is the velocity components (m/s), for all x∈[0,l].Assume thath is a constant and ζ≪h. Then (2.1) leads to (2.2)∂ζ∂t+h∂u∂x=̇0,∂u∂t+g∂ζ∂x=0.
We will consider the equation in dimensionless problem by letting U=u/gh, X=x/l, Z=ζ/h, and T=tgh/l. Substituting them into (2.2) leads to (2.3)∂Z∂T+∂U∂X=0,∂U∂T+∂Z∂X=0.
We now introduce a damping term into (2.3) to represent frictional forces due to the drag of sides of the stream, thus (2.4)∂Z∂T+∂U∂X=0,∂U∂T+∂Z∂X=-U,
with the initial conditions at t=0 and 0≤X≤1 being specified: Z=0 and U=0. The boundary conditions for t>0 are specified: Z=eit at X=0 and ∂Z/∂X=0 at X=1. The system of (2.4) is called the damped equation. We solve the damped equation by using the finite difference method. In order to solve (2.4) in [0,1]×[0,T], for convenient using u,d for U and Z, respectively, (2.5)∂u∂t+∂d∂x=-u,∂d∂t+∂u∂x=0,
with the initial conditions u=0, d=0 at t=0 and the boundary conditions d(0,t)=f(t) and ∂d/∂x=0 at x=1.


## 2.2. Dispersion Model

In a stream water quality model, the governing equations are the dynamic one-dimensional advection-dispersion-reaction equations (ADREs). A simplified representation by averaging the equation over the depth is shown in [3–5, 7, 14] as (2.6)∂C∂t+u∂C∂x=D∂2C∂x2-KC,
where C(x,t) is the concentration averaged in depth at the point x and at time t, D is the diffusion coefficient, K is the mass decay rate, and u(x,t) is the velocity component, for all x∈[0,1]. We will consider the model with following conditions. The initial condition C(x,0)=0 at t=0 for all x>0. The boundary conditions C(0,t)=C0 at x=0 and ∂C/∂x=0 at x=1, where C0 is a constant.


## 3. Numerical Experiment

The hydrodynamic model provides the velocity field and elevation of the water. Then the calculated results of the model will be input into the dispersion model which provides the pollutant concentration field.


### 3.1. Crank-Nicolson Method for the Hydrodynamic Model

We will follow the numerical techniques of [13]. To find the water velocity and water elevation from (2.5), we make the following change of variable, v=etu, and substituting them into (2.5), we have (3.1)∂v∂t+e-t∂d∂x=0,∂d∂t+e-t∂v∂x=0.
Equation (3.1) can be written in the matrix form (3.2)(vd)t+[0ete-t0](vd)x=(00).
That is, (3.3)Ut+AUx=0̅,
where (3.4)A=[0ete-t0],U=(vd),(vd)t=(∂v∂t∂d∂t),
with the initial condition d=v=0 at t=0. The left boundary conditions for x=0, t>0 are specified: d(0,t)=sint and ∂v/∂x=-etcost, and the right boundary conditions for x=1, t>0 are specified: ∂d/∂x=0 and v(0,t)=0.We now discretize (3.3) by dividing the interval [0,1] into M subintervals such that MΔx=1 and the interval [0,T] into N subintervals such that NΔt=T. We can then approximate d(xi,tn) by din, value of the difference approximation of d(x,t) at point x=iΔx and t=nΔt, where 0≤i≤M and 0≤n≤N, and similarly defined for vin and Uin. The grid points (xn,tn) are defined by xi=iΔx for all i=0,1,2,…,M and tn=nΔt for all n=0,1,2,…,N in which M and N are positive integers. Using the Crank-Nicolson method [17] to (3.3), the following finite difference equation can be obtained: (3.5)[I-14λA(Δx+∇x)]Uin+1=[I+14λA(Δx+∇x)]Uin,
where (3.6)Uin=(vindin),ΔxUin=Ui+1n-Uin,∇xUin=Uin-Ui-1n,I is the unit matrix of order 2 and λ=Δt/Δx. Applying the initial and boundary conditions given for (3.1), the general form can be obtained (3.7)Gn+1U̅n+1=EnU̅n+Fn,
where (3.8)Gn+1=[100-λ4a1n+100λ4a2n+11-λ4a2n+10000λ4a1n+1100-λ4a1n+1λ4a2n+1001-λ4a2n+10⋱⋱⋱⋱⋱⋱000λ4a1n+11-λ4a1n+100λ4a2n+1001],En=[100-λ4a1n00-λ4a2n1λ4a2n0000-λ4a1n100λ4a1n-λ4a2n001λ4a2n0⋱⋱⋱⋱⋱⋱000-λ4a1n1λ4a1n00-λ4a2n001],U̅n=(U1n+1U2n+1⋮UM-1n+1),Fn=(-λ4a1n+1sin⁡(tn+1)-λ4a1nsin⁡(tn)-λ4a2n+1Δxe-tn+1cos⁡⁡(tn+1)-λ4a2nΔxe-tncos⁡⁡(tn)00⋮00),
where a1n=etn, a2n=e-tn, and tn=nΔt for all n=0,1,2,…,N. The Crank-Nicolson scheme is unconditionally stable [15, 17].


### 3.2. The Explicit Schemes for the Advection-Diffusion-Reaction Equation

We can then approximateC(xi,tn) by Cin, the value of the difference approximation of C(x,t) at point x=iΔx and t=nΔt, where 0≤i≤M and 0≤n≤N. The grid points (xn,tn) are defined by xi=iΔx for all i=0,1,2,…,M and tn=nΔt for all n=0,1,2,…,N in which M and N are positive integers.


#### 3.2.1. Forward Time Central Space Explicit Finite Difference Scheme

Taking the forward time central space technique [17] into (2.6), we get the following discretization: (3.9)C≅Cin,∂C∂t≅Cin+1-CinΔt,∂C∂x≅Ci+1n-Ci-1n2Δx,∂2C∂x2≅Ci+1n-2Cin+Ci-1n(Δx)2,u≅Ûin.
Substituting (3.9) into (2.6), we get (3.10)Cin+1-CinΔt+Ûin(Ci+1n-Ci-1n2Δx)=D(Ci+1n-2Cin+Ci-1n(Δx)2)-KCin,
for 1≤i≤M and 0≤n≤N-1.  Leting  λ=DΔt/(Δx)2 and γin+1=(Δt/Δx)Ûin+1, (3.10) becomes (3.11)Cin+1=(12γin+λ)Ci-1n+(1-2λ-KΔt)Cin+(λ-12γin)Ci+1n.
For i=1, plugging the known value of the left boundary C0n=C0 to (3.11) in the right-hand side, we obtain(3.12)C1n+1=(12γ1n+λ)C0+(1-2λ-KΔt)C1n+(λ-12γ1n)C2n.
For i=M, substituting the approximate unknown value of the right boundary by [4], we can let CM+1n=2CMn-CM-1n and by rearranging, we obtain (3.13)CMn+1=(2λ)CM-1n+(1-4λ-KΔt-γmn)CMn.
The forward time central space scheme is conditionally stable subject to constraints in (3.10). The stability requirements for the scheme are [4, 18] (3.14)λ=DΔt(Δx)2<12,γin=UinΔtΔx<1,
where λ is the diffusion number (dimensionless) and γin is the advection number or Courant number (dimensionless). It can be obtained that the strictly stability requirements are the main disadvantage of this scheme.The finite difference formula (3.11) has been derived in [19] that the truncation error for this method is O{(Δx2),Δt}.


#### 3.2.2. Saulyev Explicit Finite Difference Scheme

The Saulyev scheme is unconditionally stable [6]. It is clear that the nonstrictly stability requirement of Saulyev scheme is the main of advantage and economical to use. Taking Saulyev technique [6] into (2.6), it can be obtained the following discretization: (3.15)C≅Cin,∂C∂t≅Cin+1-CinΔt,∂C∂x≅Ci+1n-Ci-1n+12Δx,∂2C∂x2≅Ci+1n-Cin-Cin+1+Ci-1n+1(Δx)2,u≅Ûin.
Substituting (3.15) into (2.6), we get (3.16)Cin+1-CinΔt+Ûin(Ci+1n-Ci-1n+12Δx)=D(Ci+1n-Cin-Cin+1+Ci-1n+1(Δx)2)-kCin,
for 1≤i≤M and 0≤n≤N-1. Leting λ=DΔt/(Δx)2 and γin+1=(Δt/Δx)Ûin+1, (3.16) becomes (3.17)Cin+1=(11+λ)((12γin+λ)Ci-1n+1+(1-λ-KΔt)Cin+(λ-12γin)Ci+1n).
For i=1, plugging the known value of the left boundary C0n+1=C0 to (3.17) in the right-hand side, we obtain(3.18)C1n+1=(11+λ)((12γ1n+λ)C0+(1-λ-KΔt)C1n+(λ-12γ1n)C2n).
For i=M, substituting the approximate unknown value of the right boundary by [20], we can let CM+1n=CMn+CMn+1-CM-1n+1 and by rearranging, we obtain (3.19)CMn+1=(11+(1/2)γMn)(γMnCM-1n+1+(1-KΔt-12γMn)CMn).
Using Taylor series expansions on the approximation, [21] has shown that the truncation error is O{(Δx)2+(Δt)2+(Δt/Δx)2} or O{2,2,(1/1)2}.From (3.16) to (3.19), it can be obtained that the technique does not generate the system of linear equations. It follows that the application of the technique is economical computer implementation.


## 3.1. Crank-Nicolson Method for the Hydrodynamic Model

We will follow the numerical techniques of [13]. To find the water velocity and water elevation from (2.5), we make the following change of variable, v=etu, and substituting them into (2.5), we have (3.1)∂v∂t+e-t∂d∂x=0,∂d∂t+e-t∂v∂x=0.
Equation (3.1) can be written in the matrix form (3.2)(vd)t+[0ete-t0](vd)x=(00).
That is, (3.3)Ut+AUx=0̅,
where (3.4)A=[0ete-t0],U=(vd),(vd)t=(∂v∂t∂d∂t),
with the initial condition d=v=0 at t=0. The left boundary conditions for x=0, t>0 are specified: d(0,t)=sint and ∂v/∂x=-etcost, and the right boundary conditions for x=1, t>0 are specified: ∂d/∂x=0 and v(0,t)=0.We now discretize (3.3) by dividing the interval [0,1] into M subintervals such that MΔx=1 and the interval [0,T] into N subintervals such that NΔt=T. We can then approximate d(xi,tn) by din, value of the difference approximation of d(x,t) at point x=iΔx and t=nΔt, where 0≤i≤M and 0≤n≤N, and similarly defined for vin and Uin. The grid points (xn,tn) are defined by xi=iΔx for all i=0,1,2,…,M and tn=nΔt for all n=0,1,2,…,N in which M and N are positive integers. Using the Crank-Nicolson method [17] to (3.3), the following finite difference equation can be obtained: (3.5)[I-14λA(Δx+∇x)]Uin+1=[I+14λA(Δx+∇x)]Uin,
where (3.6)Uin=(vindin),ΔxUin=Ui+1n-Uin,∇xUin=Uin-Ui-1n,I is the unit matrix of order 2 and λ=Δt/Δx. Applying the initial and boundary conditions given for (3.1), the general form can be obtained (3.7)Gn+1U̅n+1=EnU̅n+Fn,
where (3.8)Gn+1=[100-λ4a1n+100λ4a2n+11-λ4a2n+10000λ4a1n+1100-λ4a1n+1λ4a2n+1001-λ4a2n+10⋱⋱⋱⋱⋱⋱000λ4a1n+11-λ4a1n+100λ4a2n+1001],En=[100-λ4a1n00-λ4a2n1λ4a2n0000-λ4a1n100λ4a1n-λ4a2n001λ4a2n0⋱⋱⋱⋱⋱⋱000-λ4a1n1λ4a1n00-λ4a2n001],U̅n=(U1n+1U2n+1⋮UM-1n+1),Fn=(-λ4a1n+1sin⁡(tn+1)-λ4a1nsin⁡(tn)-λ4a2n+1Δxe-tn+1cos⁡⁡(tn+1)-λ4a2nΔxe-tncos⁡⁡(tn)00⋮00),
where a1n=etn, a2n=e-tn, and tn=nΔt for all n=0,1,2,…,N. The Crank-Nicolson scheme is unconditionally stable [15, 17].


## 3.2. The Explicit Schemes for the Advection-Diffusion-Reaction Equation

We can then approximateC(xi,tn) by Cin, the value of the difference approximation of C(x,t) at point x=iΔx and t=nΔt, where 0≤i≤M and 0≤n≤N. The grid points (xn,tn) are defined by xi=iΔx for all i=0,1,2,…,M and tn=nΔt for all n=0,1,2,…,N in which M and N are positive integers.


### 3.2.1. Forward Time Central Space Explicit Finite Difference Scheme

Taking the forward time central space technique [17] into (2.6), we get the following discretization: (3.9)C≅Cin,∂C∂t≅Cin+1-CinΔt,∂C∂x≅Ci+1n-Ci-1n2Δx,∂2C∂x2≅Ci+1n-2Cin+Ci-1n(Δx)2,u≅Ûin.
Substituting (3.9) into (2.6), we get (3.10)Cin+1-CinΔt+Ûin(Ci+1n-Ci-1n2Δx)=D(Ci+1n-2Cin+Ci-1n(Δx)2)-KCin,
for 1≤i≤M and 0≤n≤N-1.  Leting  λ=DΔt/(Δx)2 and γin+1=(Δt/Δx)Ûin+1, (3.10) becomes (3.11)Cin+1=(12γin+λ)Ci-1n+(1-2λ-KΔt)Cin+(λ-12γin)Ci+1n.
For i=1, plugging the known value of the left boundary C0n=C0 to (3.11) in the right-hand side, we obtain(3.12)C1n+1=(12γ1n+λ)C0+(1-2λ-KΔt)C1n+(λ-12γ1n)C2n.
For i=M, substituting the approximate unknown value of the right boundary by [4], we can let CM+1n=2CMn-CM-1n and by rearranging, we obtain (3.13)CMn+1=(2λ)CM-1n+(1-4λ-KΔt-γmn)CMn.
The forward time central space scheme is conditionally stable subject to constraints in (3.10). The stability requirements for the scheme are [4, 18] (3.14)λ=DΔt(Δx)2<12,γin=UinΔtΔx<1,
where λ is the diffusion number (dimensionless) and γin is the advection number or Courant number (dimensionless). It can be obtained that the strictly stability requirements are the main disadvantage of this scheme.The finite difference formula (3.11) has been derived in [19] that the truncation error for this method is O{(Δx2),Δt}.


### 3.2.2. Saulyev Explicit Finite Difference Scheme

The Saulyev scheme is unconditionally stable [6]. It is clear that the nonstrictly stability requirement of Saulyev scheme is the main of advantage and economical to use. Taking Saulyev technique [6] into (2.6), it can be obtained the following discretization: (3.15)C≅Cin,∂C∂t≅Cin+1-CinΔt,∂C∂x≅Ci+1n-Ci-1n+12Δx,∂2C∂x2≅Ci+1n-Cin-Cin+1+Ci-1n+1(Δx)2,u≅Ûin.
Substituting (3.15) into (2.6), we get (3.16)Cin+1-CinΔt+Ûin(Ci+1n-Ci-1n+12Δx)=D(Ci+1n-Cin-Cin+1+Ci-1n+1(Δx)2)-kCin,
for 1≤i≤M and 0≤n≤N-1. Leting λ=DΔt/(Δx)2 and γin+1=(Δt/Δx)Ûin+1, (3.16) becomes (3.17)Cin+1=(11+λ)((12γin+λ)Ci-1n+1+(1-λ-KΔt)Cin+(λ-12γin)Ci+1n).
For i=1, plugging the known value of the left boundary C0n+1=C0 to (3.17) in the right-hand side, we obtain(3.18)C1n+1=(11+λ)((12γ1n+λ)C0+(1-λ-KΔt)C1n+(λ-12γ1n)C2n).
For i=M, substituting the approximate unknown value of the right boundary by [20], we can let CM+1n=CMn+CMn+1-CM-1n+1 and by rearranging, we obtain (3.19)CMn+1=(11+(1/2)γMn)(γMnCM-1n+1+(1-KΔt-12γMn)CMn).
Using Taylor series expansions on the approximation, [21] has shown that the truncation error is O{(Δx)2+(Δt)2+(Δt/Δx)2} or O{2,2,(1/1)2}.From (3.16) to (3.19), it can be obtained that the technique does not generate the system of linear equations. It follows that the application of the technique is economical computer implementation.


## 3.2.1. Forward Time Central Space Explicit Finite Difference Scheme

Taking the forward time central space technique [17] into (2.6), we get the following discretization: (3.9)C≅Cin,∂C∂t≅Cin+1-CinΔt,∂C∂x≅Ci+1n-Ci-1n2Δx,∂2C∂x2≅Ci+1n-2Cin+Ci-1n(Δx)2,u≅Ûin.
Substituting (3.9) into (2.6), we get (3.10)Cin+1-CinΔt+Ûin(Ci+1n-Ci-1n2Δx)=D(Ci+1n-2Cin+Ci-1n(Δx)2)-KCin,
for 1≤i≤M and 0≤n≤N-1.  Leting  λ=DΔt/(Δx)2 and γin+1=(Δt/Δx)Ûin+1, (3.10) becomes (3.11)Cin+1=(12γin+λ)Ci-1n+(1-2λ-KΔt)Cin+(λ-12γin)Ci+1n.
For i=1, plugging the known value of the left boundary C0n=C0 to (3.11) in the right-hand side, we obtain(3.12)C1n+1=(12γ1n+λ)C0+(1-2λ-KΔt)C1n+(λ-12γ1n)C2n.
For i=M, substituting the approximate unknown value of the right boundary by [4], we can let CM+1n=2CMn-CM-1n and by rearranging, we obtain (3.13)CMn+1=(2λ)CM-1n+(1-4λ-KΔt-γmn)CMn.
The forward time central space scheme is conditionally stable subject to constraints in (3.10). The stability requirements for the scheme are [4, 18] (3.14)λ=DΔt(Δx)2<12,γin=UinΔtΔx<1,
where λ is the diffusion number (dimensionless) and γin is the advection number or Courant number (dimensionless). It can be obtained that the strictly stability requirements are the main disadvantage of this scheme.The finite difference formula (3.11) has been derived in [19] that the truncation error for this method is O{(Δx2),Δt}.


## 3.2.2. Saulyev Explicit Finite Difference Scheme

The Saulyev scheme is unconditionally stable [6]. It is clear that the nonstrictly stability requirement of Saulyev scheme is the main of advantage and economical to use. Taking Saulyev technique [6] into (2.6), it can be obtained the following discretization: (3.15)C≅Cin,∂C∂t≅Cin+1-CinΔt,∂C∂x≅Ci+1n-Ci-1n+12Δx,∂2C∂x2≅Ci+1n-Cin-Cin+1+Ci-1n+1(Δx)2,u≅Ûin.
Substituting (3.15) into (2.6), we get (3.16)Cin+1-CinΔt+Ûin(Ci+1n-Ci-1n+12Δx)=D(Ci+1n-Cin-Cin+1+Ci-1n+1(Δx)2)-kCin,
for 1≤i≤M and 0≤n≤N-1. Leting λ=DΔt/(Δx)2 and γin+1=(Δt/Δx)Ûin+1, (3.16) becomes (3.17)Cin+1=(11+λ)((12γin+λ)Ci-1n+1+(1-λ-KΔt)Cin+(λ-12γin)Ci+1n).
For i=1, plugging the known value of the left boundary C0n+1=C0 to (3.17) in the right-hand side, we obtain(3.18)C1n+1=(11+λ)((12γ1n+λ)C0+(1-λ-KΔt)C1n+(λ-12γ1n)C2n).
For i=M, substituting the approximate unknown value of the right boundary by [20], we can let CM+1n=CMn+CMn+1-CM-1n+1 and by rearranging, we obtain (3.19)CMn+1=(11+(1/2)γMn)(γMnCM-1n+1+(1-KΔt-12γMn)CMn).
Using Taylor series expansions on the approximation, [21] has shown that the truncation error is O{(Δx)2+(Δt)2+(Δt/Δx)2} or O{2,2,(1/1)2}.From (3.16) to (3.19), it can be obtained that the technique does not generate the system of linear equations. It follows that the application of the technique is economical computer implementation.


## 4. The Accuracy of the Hydrodynamic Approximation

It is not hard to find the analytical solutiond(x,t) in (2.5) with f(t)=sint. By changing of variables, d(x,t)=eit𝒟(x) and u(x,t)=eit𝒰(x) for some 𝒟(x),𝒰(x)∈C02[0,1] and substitute into (2.5). Using a separable variables technique, we can obtain d(1,t) a solution [14] (4.1)d(1,t)=sin⁡tcos⁡⁡βcosh⁡⁡α-cos⁡⁡tsin⁡βsinh⁡⁡αcos⁡2⁡βcosh⁡2⁡α+sin2⁡βsinh⁡2⁡α,
where α=21/4cos(3π/8) and β=21/4sin(3π/8). Anyhow, it is not easy to find the analytical solution u(x,t) of (2.5). We use the solution d(1,t) obtained in (4.1) to verify to our approximate solution obtained by the Crank-Nicolson method (3.7). Actually when using the Crank-Nicolson method, we get the approximate solution both d(x,t) and u(x,t). We assume that when we get a good approximation for d(x,t) this implied that the method gives a good approximation for u(x,t). The verification of the approximate solution d̅(1,t) is shown in Table 1 and Figure 2.Table 1
The error defined byerror(T)=max|d(1,t)-d̅(1,t)| for all T-π≤t≤T+π at fixed Δt/Δx=0.25.
TΔxΔtd(1,t)d̅(1,t)error(T)200.2000.050001.501581.456590.044990.1000.025000.093920.388560.294650.0500.01250−0.25831−0.402440.144130.0250.00625−0.33323−0.404330.07110300.2000.050001.496661.451650.045010.1000.02500−0.09358−0.388210.294630.0500.012500.257980.401980.144010.0250.006250.332900.403890.07099400.2000.050001.501461.456440.045020.1000.025000.093240.387870.294630.0500.01250−0.25765−0.401650.144010.0250.00625−0.33257−0.403560.07099Figure 2
Comparison of analytical solution for height of water elevation with results obtained by numerical technique at the end point of the domain.Figure2 shows the comparison between the analytical solutions d(1,t) and the approximate solutions d̅(1,t) only at the end of the domain. Table 1 shows that an estimate of the maximum error is less than 7.0%.Unfortunately, the analytical solutions of hydrodynamic model could not found over entire domain. This implies that the analytical solutions of dispersion model could not carry out at any points on the domain as well.


## 5. Application to the Stream Water Quality Assessment Problem

Suppose that the measurement of pollutant concentrationC (Kg/m3) in a uniform stream at time t (sec) is considered. A stream is aligned with longitudinal distance, 1.0 km total length, and 1.0 m depth. There is a plant which discharges waste water into the stream, and the pollutant concentration at discharge point is C(0,t)=C0=1 Kg/m3 at x=0 for all t>0 and C(x,0)=0 Kg/m3 at t=0. The elevation of water at the discharge point can be described as a function d(0,t)=f(t)=sint m for all t>0, and the elevation is not changed at x=1.0 km. The physical parameters of the stream system are diffusion coefficient D=0.05 m2/s, and a first-order reaction rate 10-5 s−1. In the analysis conducted in this study, meshes the stream into 20, 40, 80, and 160 elements with Δx=0.05,0.025,0.0125,0.00625, respectively, and time increment are 3.2, 1.6, 0.8, 0.4, 0.2, 0.1 s with Δt=0.01,0.005,0.0025,0.00125,0.000625,0.0003125, respectively. Using (3.7), it can be obtained the water velocity u(x,t) on Table 2 and Figure 3. Next, the approximate water velocity can be plugged into the finite difference equations of FTCS and Saulyev schemes on (3.11)–(3.13) and (3.17)–(3.19), respectively. The comparison of concentration of the FTCS method and Saulyev method is presented in Figure 5 for two different instants. We then have the stabilities of both schemes for each Δx and Δt in Table 5 that are consistence with (3.14). The approximation of pollutant concentration C of both schemes are shown in Tables 3 and 4 and Figure 4. The concentration along a stream at only 24 min with varied diffusion coefficients is shown in Figure 6. These imply that the Peclet number was Pe=(UinΔx/D)<10, which indicated the stream system was advection not dominated [15].Table 2
The velocity of water flowu(x,t).
tx=0x=0.1x=0.2x=0.3x=0.4x=0.5x=0.6x=0.7x=0.8x=0.9x=1.010−0.5478−0.5695−0.5666−0.5199−0.4711−0.4061−0.3352−0.2571−0.1747−0.08780.0000201.31011.22131.12321.01600.89310.75040.60970.45950.30990.15540.000030−0.4468−0.3731−0.3085−0.2527−0.2057−0.1651−0.1252−0.0875−0.0573−0.02760.000040−1.0361−0.9898−0.9258−0.8459−0.7513−0.6439−0.5258−0.4004−0.2711−0.13670.0000501.09390.99180.88670.77910.67000.55940.44790.33560.22330.11140.0000Table 3
The pollutant concentrationC(x,t) of FTCS scheme, Δx=0.0125 (40 m), Δt=0.0003125 (0.1 s).
tx=0x=0.1x=0.2x=0.3x=0.4x=0.5x=0.6x=0.7x=0.8x=0.9x=1.0101.00000.32310.11540.05170.02920.01940.01410.01090.00900.00790.0075201.00000.99960.99710.98730.95750.88730.75880.57740.38280.23240.1720301.00000.99810.99560.99100.98300.97060.95330.93270.91220.89660.8903401.00000.92600.91660.91500.91450.91410.91380.91360.91350.91340.9133501.00000.99930.99660.99020.97880.96330.94680.93270.92310.91790.9162Table 4
The pollutant concentrationC(x,t) of Saulyev scheme, Δx = 0.0125 (40 m), Δt=0.005 (1.6 s).
tx=0x=0.1x=0.2x=0.3x=0.4x=0.5x=0.6x=0.7x=0.8x=0.9x=1.0101.00000.32970.12250.05700.03270.02120.01440.00960.00590.00280.0000201.00000.99950.99700.98780.96190.90200.79080.62510.42380.21360.0000301.00000.99720.99190.97980.95360.90260.81360.67600.48590.25300.0000401.00000.31280.20310.16240.13210.10520.08070.05810.03720.01790.0000501.00000.99130.96310.89910.78790.63420.46100.29880.16940.07530.0000Table 5
The stability of FTCS and Saulyev schemes,γ=max{γin:1≤i≤M,0≤n≤N-1}.
ΔxΔtλγFTCS schemeSaulyev scheme0.01250.0100003.21.0665UnstableStable0.0050001.60.5332UnstableStable0.0025000.80.2667UnstableStable0.0012500.40.1333StableStable0.0006250.20.0080StableStable0.01250.0100003.21.0665UnstableStable0.0050001.60.5332UnstableStable0.0025000.80.2667UnstableStable0.0012500.40.1333StableStable0.0006250.20.0080StableStable0.01250.0100003.21.0665UnstableStable0.0050001.60.5332UnstableStable0.0025000.80.2667UnstableStable0.0012500.40.1333StableStable0.0006250.20.0080StableStableFigure 3
The water velocityu(x,t).Figure 4
The pollutant concentrationC(x,t).The comparison of concentration at two different time instants of the FTCS and Saulyev methods.
(a)(b)The pollutant concentrationC(x,t) for each diffusion coefficients at t=24 min.
(a)(b)(c)(d)


## 6. Discussion and Conclusions

The approximation of the pollutant concentrations from the FTCS technique is shown in Tables3 and 5; it can be concluded that stability requirements is one of the disadvantages of the technique. The real-world problems require a small amount of time interval in obtaining accurate solutions. We can see that the FTCS scheme is not good agreement for real application. In Table 4, it can be obtained that the Saulyev technique has an advantage over compared to FTCS. It is unconditionally stable, easy, and economical to implement.By Figure6, we can see that the diffusion coefficients of pollutant matter can reduce the concentration in a nonuniform stream. If sewage effluent with a low diffusion coefficient has discharged into a nonuniform flow stream, then the water quality will be lower than a discharging of high diffusion coefficients of other pollutant matters.In this paper, it can be combined the hydrodynamic model and the convection-diffusion-reaction equation to approximate the pollutant concentration in a stream when the current which reflects water in the stream is not uniform. The technique developed in this paper the response of the stream to the two different external inputs: the elevation of water and the pollutant concentration at the discharged point. The Saulyev technique can be used in the dispersion model since the scheme is very simple to implement. By the Saulyev finite difference formulation, we obtain that the proposed technique is applicable and economical to be used in the real-world problem as aresult of the simplicity of programming and the straight forwardness of the implementation. It is also possible to find tentative better locations and the periods of time of the different discharged points to a stream.



---

*Source: 491317-2011-09-06.xml*