Modeling pollution transmission in rivers is an important subject in environmental engineering studies. Numerical approaches to modeling pollution transmission in rivers are useful tools for managing the water quality. The advection-dispersion equation is the governing equation in the transport of pollution in rivers. Recently, due to advances in fractional calculus in engineering modeling, the simulation of pollution transmission in rivers has been improved using the fractional derivative approach. In this study, by solving the fractional advection-dispersion equation (FRADE), a numerical model was developed for the simulation of pollution transmission in rivers with stagnant zones. To this purpose, both terms of the FRADE equation (advection and fractional dispersion) were discretized separately and the results were connected using a time-splitting technique. The analytical solution of a modified advection-dispersion equation (MADE) model and observed data from the Severn River in the UK were used to demonstrate the model capabilities. Results indicated that there is a good agreement between observed data, the analytical solution of the MADE model, and the results of the developed numerical model. The developed numerical model can accurately simulate the long-tailed dispersion processes in a natural river.

## NOMENCLATURE

• dead zone effect

•
• ability to self-clean

•

•
• C

concentration

•
• DL

longitudinal dispersion coefficient

•

•

•
• TSM

transient storage zone

## INTRODUCTION

1
where is the longitudinal dispersion coefficient, C is the concentration of pollution along the river and u is the mean velocity of river flow. The analytical solution of this equation leads to a Gaussian distribution but observed data from the field studies have shown that the concentration of tracer along the river does not follow the theoretical solution (Pujol & Sanchez-Cabeza 1999). Geometrical properties of rivers such as riffles, pools and dead zones, as shown in Figure 1, influence the concentration profile of tracer in the river and cause its shape to change to that shown in Figure 2. Figure 2 shows that the stagnant zones cause a long tail in the concentration profile of tracers. This means that some of the mass of the tracer is temporarily held in the dead zones and as time passes will return to the river flow (Mirbagheri et al. 2009).
Figure 1

Stagnant zone conditions along the river path (a) plan view, (b) side view (Singh 2003).

Figure 1

Stagnant zone conditions along the river path (a) plan view, (b) side view (Singh 2003).

Figure 2

Non-Gaussian profile of the tracer concentration (Singh 2008).

Figure 2

Non-Gaussian profile of the tracer concentration (Singh 2008).

By observing non-uniformity in the concentration of tracer in rivers, investigators have tried to improve the accuracy of pollution transmission modeling using the modification of the ADE. Singh (2003) proposed the modified advection-dispersion equation (MADE) method to improve the precision of modeling pollution transmission in rivers. The formulation of MADE is presented in Equation (2).
2
in which is the dead zone effect and is the ability of rivers to self-clean. The other parameters have been defined previously. The transient storage zone (TSM) model has been widely used for modeling pollution transmission in rivers with dead zone conditions (Singh 2003; Deng et al. 2004; Singh 2008). Seo & Cheong (2001), using theoretical approaches, estimated the required parameters of TSM and found that the moment-based approaches are the most effective to estimate the value of TSM parameters. Deng et al. (2004) proposed the fractional advection-dispersion equation (FRADE) model to simulate pollution transmission in a river with dead zone conditions. They applied the FRADE model to simulate pollution transmission in rivers and found that the F factor has a value between 1.4 and 2. The formulation of the FRADE model is given in Equation (3):
3

## MATERIALS AND METHODS

In this study, a computational code has been written corresponding to the FRADE model, a subroutine has been written for the analytical solution of MADE, and the performance of the FRADE model is evaluated.

### Analytical solution of MADE

Based on recommendations of Singh (2003) and assuming velocity of flow and dispersion coefficient as constant values as well as a sudden injection of tracer, the analytical solution of MADE can be derived as Equation (4):
4a
4b
4c
where ua is the apparent flow velocity and Da is the corrected dispersion coefficient.

### Proposed numerical model for FRADE

To provide the numerical model for FRADE, a time splitting method has been used. In this way, with regard to Equation (3), two terms of the FRADE model (advection and fractional dispersion) are discrete and their results are combined using a time splitting technique. To discretize the advection term, the Fromm scheme is used and the discrete form of advection is derived as Equation (5):
5
where . To discretize the dispersion term of the FRADE model, the application of a suitable numerical scheme is considered. Several discretization methods have been proposed to this end but the more accepted Grunewald's method is used. Based on Grunewald's definition, the partial differentiation equation related to the dispersion term is defined as a finite series as in Equation (6):
6
where , N is a natural number and is a gamma function. Oldham & Spanier (2006) stated that Equation (6) has a better convergence than using the Grunewald definition. To calculate the concentration of tracer in a desired place, Equation (6) should be evaluated. A simple way to achieve this is by using the Deng-Singh-Bengtsson definition (Equation (7)):
7
where n is the time and wjF = (−1)j is the weighted factor and . Values of for different types of other differentiations are given in Table 1. The highest value of is 0.064 when j = 3 and F = 1.6. The finite series (Equation (7)) can be rewritten as Equation (8).
8a
8b
In Equation (8), the first term on the right-hand side involves the Gaussian curve and FT1 is an attempt to model the fractional tail (FT). In this study, the fully implicit scheme is developed and used as Equation (9):
9
where m = 1,2,3,. . . , N − 1 and j = i+1-m. Equation (9) can be rewritten as Equation (10):
10
where and R = Cmn + .
For m=1 to N1, a three-diagonal matrix is constructed and is solved using the Thomas algorithm. To apply the proposed numerical model in computer, a script was written in Matlab software and a flowchart for calibration of the proposed numerical model is given in Figure 3.
Table 1

Changes in fractional binomial coefficients with memory length j and factor F results

Factor F
Memory length j2.0001.9001.8001.7001.6001.5001.4001.3001.2001.1001.000
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
−2.000 −1.90 −1.80 −1.70 −1.60 −1.50 −1.40 −1.30 −1.20 −1.10 −1.00
1.000 0.855 0.720 0.595 0.480 0.375 0.280 0.195 0.120 0.055 0.000
0.000 0.029 0.048 0.060 0.064 0.063 0.056 0.046 0.032 0.017 0.000
0.000 0.008 0.014 0.019 0.022 0.023 0.022 0.019 0.014 0.008 0.000
0.000 0.003 0.006 0.009 0.011 0.012 0.012 0.010 0.008 0.005 0.000
0.000 0.002 0.003 0.005 0.006 0.007 0.007 0.006 0.005 0.003 0.000
0.000 0.001 0.002 0.003 0.004 0.004 0.005 0.004 0.004 0.002 0.000
0.000 0.001 0.001 0.002 0.003 0.003 0.003 0.003 0.003 0.002 0.000
0.000 0.000 0.001 0.001 0.002 0.002 0.002 0.002 0.002 0.001 0.000
10 0.000 0.000 0.001 0.001 0.001 0.002 0.002 0.002 0.002 0.001 0.000
11 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.000
13 0.000 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.000
15 0.000 0.000 0.000 0.000 0.000 0.001 0.001 0.001 0.001 0.000 0.000
Factor F
Memory length j2.0001.9001.8001.7001.6001.5001.4001.3001.2001.1001.000
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
−2.000 −1.90 −1.80 −1.70 −1.60 −1.50 −1.40 −1.30 −1.20 −1.10 −1.00
1.000 0.855 0.720 0.595 0.480 0.375 0.280 0.195 0.120 0.055 0.000
0.000 0.029 0.048 0.060 0.064 0.063 0.056 0.046 0.032 0.017 0.000
0.000 0.008 0.014 0.019 0.022 0.023 0.022 0.019 0.014 0.008 0.000
0.000 0.003 0.006 0.009 0.011 0.012 0.012 0.010 0.008 0.005 0.000
0.000 0.002 0.003 0.005 0.006 0.007 0.007 0.006 0.005 0.003 0.000
0.000 0.001 0.002 0.003 0.004 0.004 0.005 0.004 0.004 0.002 0.000
0.000 0.001 0.001 0.002 0.003 0.003 0.003 0.003 0.003 0.002 0.000
0.000 0.000 0.001 0.001 0.002 0.002 0.002 0.002 0.002 0.001 0.000
10 0.000 0.000 0.001 0.001 0.001 0.002 0.002 0.002 0.002 0.001 0.000
11 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.000
13 0.000 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.000
15 0.000 0.000 0.000 0.000 0.000 0.001 0.001 0.001 0.001 0.000 0.000
Figure 3

Flowchart for training proposed numerical model.

Figure 3

Flowchart for training proposed numerical model.

Based on this flowchart, first the advection term is solved using the Fromm scheme and then the result is entered into the dispersion section concurrent with simulation; then the final results are derived. If the results are not acceptable, the model is trained again using a new value for dispersion, time and spatial step sizes. To validate the performance of the proposed model, a comparison is made with the analytical solution of MADE and modeling pollution transmission in the River Severn.

### Site description

Atkinson & Davis (2000) studied the mechanism of pollutant transmission in the River Severn. They considered 14 km of the river and six sampling stations measuring the concentration of tracer. They investigated the effect of hydraulic and geometric properties of the river such as bed form, reverse flow and dead zones on the profile of the tracer concentration. Figure 4 shows a schematic of the River Severn and the location of sampling stations (Davis & Atkinson 1999; Davis et al. 1999). As mentioned in a description of the flowchart, one of the inputs is . To provide a good estimation of this parameter and reduce the iteration of the trial and error process when training the model, was calculated using the routing concentration profile method proposed by Disley et al. (2015). Using this method includes a number of steps as follows: considering the initial value for , calculating the concentration profile at the downstream station using the upstream concentration profile and , and performing a comparison between the calculated profile and the measured profile. If the calculated profile does not suitably cover the measured profile, the process will be repeated until the calculated profile has a good coverage of the measured profile.
Figure 4

Schematic map of the River Severn and sampling stations (Atkinson & Davis 2000).

Figure 4

Schematic map of the River Severn and sampling stations (Atkinson & Davis 2000).

## RESULTS AND DISCUSSION

The results of running the prepared script for the Fromm scheme are given in Figure 5. As shown in this figure, after passing about 1,600 m, the concentration peak of the tracer does not decay. Evaluating the performance of the scheme shows that this model is suitable to model pure advection. For evaluating the pure advection term, , , .
Figure 5

Results of solving the pure advection term by the Fromm scheme.

Figure 5

Results of solving the pure advection term by the Fromm scheme.

The results of running the script for a pure partial differential equation considering =5 L2/s, Δx = 8 m, Δt = 2 m and α = 1.85 are shown in Figure 6. As shown in this figure, as time passes, the skewness of the curves increases.
Figure 6

Results of solving the pure dispersion term with the fully implicit scheme.

Figure 6

Results of solving the pure dispersion term with the fully implicit scheme.

The results of the MADE analytical solution are presented in Figure 7. These results were obtained considering u = 1.1 m/s and DL = 10.2 L2/s. As shown in this figure, by decreasing the concentration peak, the skewness of the curves increases.
Figure 7

Results of the MADE analytical solution.

Figure 7

Results of the MADE analytical solution.

A comparison of the results of the proposed numerical model for the FRADE model and the MADE analytical solution is given in Figure 8. For both models, Δx = 10 m, Δt = 2 s, =2.5, ɛ = 0.25. This result was taken at T = 270 s. As shown in this figure, the performance of the developed numerical model is suitable in comparison with the results of the MADE analytical solution.
Figure 8

Results of the MADE analytical solution and the developed numerical model.

Figure 8

Results of the MADE analytical solution and the developed numerical model.

To validate the performance of the developed numerical model for practical purposes, a simulation of tracer transmission in the River Severn was considered. As mentioned earlier, in order to apply the model, it first has to be trained. The dispersion coefficient method was used to reduce the trial and error iterations. As presented in this paper, during training the developed numerical model dispersion coefficient and F factor were derived as 59.2 and 1.68, respectively. Results of derived values of and the F factor for other stations are given in Table 2. As seen from Table 2, by increasing the tail part of the profile of tracer concentration, the values of fractional order are decreased. As presented in Table 2, at the primary stations where the dead zones do not have a full impact on the tracer concentration profile the fractional order obtained was equal to 1.68 and by increasing the effects of dead zones on the tracer concentration profile the order of fractional dispersion was decreased to 1.48. Performance of the developed model trained with published information at station A is shown in Figure 9 and the results of testing the model for the published information at other stations are given in Figure 10. As shown in these figures, performance of the developed model for the testing (validation) stage is suitable.
Table 2

Value of input parameters for simulation of pollution transmission in the River Severn

StationSpatial step (m)Time step (s)Fractional order (F factor)Dispersion coefficient
1.68 59.2
1.68 19.2
1.68 12.8
1.59 26.3
1.53 37.4
1.48 29.7
StationSpatial step (m)Time step (s)Fractional order (F factor)Dispersion coefficient
1.68 59.2
1.68 19.2
1.68 12.8
1.59 26.3
1.53 37.4
1.48 29.7
Figure 9

The results of calibration of numerical model for stations (A).

Figure 9

The results of calibration of numerical model for stations (A).

Figure 10

The results of validation of the numerical model.

Figure 10

The results of validation of the numerical model.

## CONCLUSION

Modeling pollution transmission in rivers is an important part of environmental engineering studies. Stagnant zones existing in rivers temporarily hold some of the pollutant, and as time passes it is returned into the main river flow. This temporary holding of pollution causes a long tail in the concentration profile along rivers. A long tail in a concentration profile can not be modeled using classical ADE as the governing equation of pollution transmission in rivers. Due to advances in fractional calculus in most areas of physics and engineering, investigators have used this technique for improving the precision of modeling engineering problems. In this study, pollution transmission in a river was modeled using the FRADE model. In this model, the dispersion term is modeled by a fractional partial differential equation. To prepare a numerical model based on the FRADE approach, the advection term was discretized using the finite volume method and the fractional dispersion term was discretized with a fully implicit central scheme. The main parameter in terms of the FRADE model is determining the order of the fractional dispersion (F) term. Results of calibration of the developed numerical model for the River Severn indicated that the order of fractional dispersion of 1.68 can accurately model the long tail of the concentration profile. By moving the tracer along the river and increasing the tail of the concentration profile, the values of the fractional order are decreased.

## REFERENCES

Altunkaynak
A.
2016
.
Journal of Hydro-Environment Research
12
,
105
116
.
doi:http://dx.doi.org/10.1016/j.jher.2016.05.001.
Atkinson
T. C.
Davis
P. M.
2000
.
Hydrology and Earth System Sciences
4
,
345
353
.
Azamathulla
H. M.
Wu
F.-C.
2011
.
Applied Soft Computing
11
,
2902
2905
.
doi:http://dx.doi.org/10.1016/j.asoc.2010.11.026.
Azamathulla
H. M.
Haghiabi
A. H.
Parsaie
A.
2016
.
Water Science and Technology: Water Supply
16
,
1002
1016
.
doi:10.2166/ws.2016.014.
Cho
Y.
2016
.
Water Science and Technology: Water Supply
16
,
703
714
.
doi:10.2166/ws.2015.180.
Davis
P. M.
Atkinson
T. C.
1999
.
Hydrology and Earth System Sciences
4
,
373
381
.
doi:10.5194/hess-4-373-2000.
Davis
P. M.
Atkinson
T. C.
Wigley
T. M. L.
1999
.
Hydrology and Earth System Sciences
4
,
355
371
.
doi:10.5194/hess-4-355-2000.
Dehdar-behbahani
S.
Parsaie
A.
2016
.
Alexandria Engineering Journal
55
,
467
473
.
doi:10.1016/j.aej.2016.01.006.
Deng
Z.-Q.
Singh
V. P.
Bengtsson
L.
2004
.
Journal of Hydraulic Engineering
130
,
422
431
.
doi:10.1061/(asce)0733-9429(2004)130:5(422).
Disley
T.
Gharabaghi
B.
Mahboubi
A. A.
McBean
E. A.
2015
.
Hydrological Processes
29
,
161
172
.
doi:10.1002/hyp.10139.
M.
O.
Seifollahi-Aghmiuni
S.
Loáiciga
H. A.
2016
.
Journal of Irrigation and Drainage Engineering
142
(
11
),
04016050
.
doi:10.1061/(asce)ir.1943-4774.0001083.
Guan
X.-j.
Liang
S.-x.
Meng
Y.
2016
.
Water Science and Technology: Water Supply
.
doi:10.2166/ws.2016.057.
Haghiabi
A. H.
2016
.
Journal of Earth System Science
125
,
985
995
.
doi:10.1007/s12040-016-0708-8.
Kashefipour
S. M.
Falconer
R. A.
2002
.
Water Research
36
,
1596
1608
.
doi:http://dx.doi.org/10.1016/S0043-1354(01)00351-7.
B.
Jafari
H.
Baleanu
D.
2013a
.
The European Physical Journal Special Topics
222
,
1805
1812
.
doi:10.1140/epjst/e2013-01965-1.
B.
Naseri
A. A.
Jafari
H.
A.
Baleanu
D.
2013b
.
Computers & Mathematics with Applications
66
,
785
794
.
doi:http://dx.doi.org/10.1016/j.camwa.2013.01.002.
Mirbagheri
S.
Abaspour
M.
Zamani
K.
2009
Mathematical modeling of water quality in river systems
,
European Water
27/28
,
31
41
.
M.
Bonakdari
H.
2016
.
Journal of Pipeline Systems Engineering and Practice
8
(
1
),
06016003
.
doi:10.1061/(asce)ps.1949-1204.0000249.
M.
Sattar
A. A.
2015
.
Water Resources Management
29
,
2205
2219
.
doi:10.1007/s11269-015-0936-8.
M.
Tafarojnoruz
A.
2016
.
Environmental Earth Sciences
75
,
1
12
.
doi:10.1007/s12665-015-4877-6.
M.
Zahiri
A.
2015
.
Journal of Hydrologic Engineering
20
,
04015035
.
doi:10.1061/(asce)he.1943-5584.0001185.
M.
A.
Lim
S. Y.
2016
.
Ocean Engineering
111
,
128
135
.
Noori
R.
Karbassi
A.
Farokhnia
A.
Dehghani
M.
2009
.
Environmental Engineering Science
26
,
1503
1510
.
doi:10.1089/ees.2008.0360.
Noori
R.
Karbassi
A.
Khakpour
A.
Shahbazbegian
M.
H. M. K.
Vesali-Naseh
M.
2011a
.
Environmental Modeling & Assessment
17
,
411
420
.
doi:10.1007/s10666-011-9302-2.
Noori
R.
Karbassi
A. R.
H.
Vesali-Naseh
M.
Sabahi
M. S.
2011b
.
Environmental Progress & Sustainable Energy
30
,
439
449
.
doi:10.1002/ep.10478.
Noori
R.
Karbassi
A. R.
A.
Han
D.
Zokaei-Ashtiani
M. H.
Farokhnia
A.
Gousheh
M. G.
2011c
.
Journal of Hydrology
401
,
177
189
.
doi:http://dx.doi.org/10.1016/j.jhydrol.2011.02.021.
Noori
R.
Deng
Z.
A.
Kachoosangi
F. T.
2015
Journal of Hydraulic Engineering
142
(
1
),
04015039
.
doi:10.1061/(asce)hy.1943-7900.0001062.
Oldham
K. B.
Spanier
J.
2006
The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order. Dover Publications
.
Parsaie
A.
Haghiabi
A. H.
2015a
.
International Journal of Scientific Research in Environmental Sciences
3
,
199
207
.
doi:10.12983/ijsres-2015-p0199-0207
.
Parsaie
A.
Haghiabi
A. H.
2015b
.
Applied Water Science
1
10
.
doi:10.1007/s13201-015-0319-6.
Parsaie
A.
Haghiabi
A. H.
2015c
.
Modeling Earth Systems and Environment
1
,
1
8
.
Parsaie
A.
Haghiabi
A. H.
2016
.
Sustainable Water Resources Management
1
8
.
doi:10.1007/s40899-016-0055-6.
Parsaie
A.
Najafian
S.
Shamsi
Z.
2016
.
Modeling Earth Systems and Environment
2
,
1
9
.
doi:10.1007/s40808-016-0207-6.
Pujol
L.
Sanchez-Cabeza
J. A.
1999
.
Journal of Environmental Radioactivity
45
,
39
57
.
doi:http://dx.doi.org/10.1016/S0265-931X(98)00075-7.
Qishlaqi
A.
Kordian
S.
Parsaie
A.
2016
.
Applied Water Science
1
6
.
doi:10.1007/s13201-016-0409-0.
Rathnayake
U.
2015a
Enhanced water quality modelling for optimal control of drainage systems under SWMM constraint handling approach
.
Asian Journal of Water, Environment and Pollution
12
,
81
85
.
Rathnayake
U.
2015b
Multi-objective optimization of combined sewer systems using SWMM 5.0
.
Journal of Civil Engineering and Architecture Research
2
,
985
993
.
Rathnayake
U. S.
Tanyimboh
T. T.
2015
.
Water Resources Management
29
,
2715
2731
.
doi:10.1007/s11269-015-0965-3.
H.
S. A.
E.
M. M.
2009
.
Expert Systems with Applications
36
,
8589
8596
.
Sahin
S.
2014
.
Environmental Processes
1
,
277
285
.
doi:10.1007/s40710-014-0018-6.
Sattar
A. M. A.
Gharabaghi
B.
2015
.
Journal of Hydrology
524
,
587
596
.
doi:http://dx.doi.org/10.1016/j.jhydrol.2015.03.016.
Seo
I. W.
Cheong
T. S.
1998
.
Journal of Hydraulic Engineering
124
,
25
32
.
doi:10.1061/(asce)0733-9429(1998)124:1(25).
Seo
I. W.
Cheong
T. S.
2001
.
Journal of Hydraulic Engineering
127
,
453
465
.
doi:10.1061/(asce)0733-9429(2001)127:6(453).
Shen
C.
Niu
J.
Anderson
E. J.
Phanikumar
M. S.
2010
.
Advances in Water Resources
33
,
615
623
.
Singh
S. K.
2003
.
Journal of Hydraulic Engineering
129
,
470
473
.
doi:10.1061/(asce)0733-9429(2003)129:6(470).
Singh
S. K.
2008
.
Journal of Irrigation and Drainage Engineering
134
,
853
856
.
doi:10.1061/(asce)0733-9437(2008)134:6(853).
Toprak
Z. F.
Şen
Z.
Savci
M. E.
2004
.
Water Research
38
,
3139
3143
.
doi:http://dx.doi.org/10.1016/j.watres.2003.08.004.
Toprak
Z. F.
Hamidi
N.
Kisi
O.
Gerger
R.
2014
.
KSCE Journal of Civil Engineering
18
,
718
730
.
doi:10.1007/s12205-014-0089-y.
Tutmez
B.
Yuceer
M.
2013
.
Water Resources Management
27
,
3307
3318
.
doi:10.1007/s11269-013-0348-6.
Vries
D.
van den Akker
B.
Vonk
E.
de Jong
W.
van Summeren
J.
2016
.
Water Science and Technology: Water Supply
.
doi:10.2166/ws.2016.062.
Wang
D.
2016
.
Water Science and Technology: Water Supply
16
,
746
755
.
doi:10.2166/ws.2015.186.
Wang
Y.
Huai
W.
2016
.
Journal of Hydraulic Engineering
142
(
11
),
04016048
.
doi:10.1061/(asce)hy.1943-7900.0001196.
West
J. R.
Mangat
J. S.
1986
.
Estuarine, Coastal and Shelf Science
22
,
161
181
.
doi:http://dx.doi.org/10.1016/0272-7714(86)90111-3.
Xu
E.
Zhang
H.
Dong
G.
Kang
L.
Zhen
X.
2016
.
Water Science and Technology: Water Supply
16
,
817
827
.
doi:10.2166/ws.2016.005.
Zeng
Y.
Huai
W.
2014
.
Journal of Hydro-Environment Research
8
,
2
8
.
doi:http://dx.doi.org/10.1016/j.jher.2013.02.005.