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

## INTRODUCTION

*et al.*2016; Wang 2016). Modeling contaminant transmission in rivers is one aspect of studying river water quality (Rathnayake 2015a; Parsaie & Haghiabi 2015b; Guan

*et al.*2016; Xu

*et al.*2016). Studies on pollution transmission have been conducted using field measurements and numerical simulations (Parsaie & Haghiabi 2015a; Rathnayake & Tanyimboh 2015). In field studies, to define the mechanism of pollution transmission, a tracer is injected and its concentration is routed along the river (Toprak

*et al.*2004; Noori

*et al.*2011a, 2011c; Tutmez & Yuceer 2013; Toprak

*et al.*2014; Zeng & Huai 2014). In other words, to characterize the effects of hydraulic and geometric properties of rivers such as dead or stagnant zones, a profile of the concentration of the tracer at a number of sampling stations is evaluated (Rathnayake 2015b; Farhadian

*et al.*2016; Wang & Huai 2016). One of the main parameters related to studies of river pollution is longitudinal dispersion coefficient () (West & Mangat 1986; Seo & Cheong 1998, 2001; Kashefipour & Falconer 2002). This parameter is usually derived by routing the concentration of tracer along the river. Recently, advanced numerical approaches (Azamathulla

*et al.*2016; Dehdar-behbahani & Parsaie 2016; Parsaie & Haghiabi 2016; Parsaie

*et al.*2016) such as image processing (Shen

*et al.*2010; Altunkaynak 2016) and soft computing techniques such as artificial neural networks (Parsaie & Haghiabi 2015c; Vries

*et al.*2016), adaptive neuro-fuzzy inference systems (Riahi-Madvar

*et al.*2009; Najafzadeh

*et al.*2016), support vector machines (Noori

*et al.*2009; Azamathulla & Wu 2011; Noori

*et al.*2015), genetic programming (Sattar & Gharabaghi 2015), the Group Method of Data Handling (GMDH) (Najafzadeh & Sattar 2015; Najafzadeh & Zahiri 2015; Najafzadeh & Bonakdari 2016; Najafzadeh & Tafarojnoruz 2016), and multivariate adaptive regression splines (Haghiabi 2016) have been applied to estimate the in rivers. When pollution is first injected in the river the dispersion process is three-dimensional, however, one may assume it becomes one-dimensional after spreading across the river and moving downstream (Sahin 2014). The governing equation is the advection-dispersion equation (ADE). A one-dimensional ADE is presented in Equation (1) (Noori

*et al.*2011b). 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).

*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): where

*F*is the factor of the effect of dead zones. The

*F*factor can have a value between 1 and 2. If

*F*has a value equal to 2, the FRADE model returns to the classical ADE. Singh (2003) applied the TSM, MADE and FRADE models to simulate pollution transmission in rivers and made the following observations. First, the MADE model needs few parameters and using this model is simple and, also, this model has an analytical solution which can easily be applied. Second, estimating the TSM model's parameters needs the application of theoretical and optimization approaches, and then using its analytical solution is practical. Using the FRADE model requires the application of a numerical solution and evaluating the

*F*factor. This model needs the use of optimization methods to determine the optimal value of

*F*. By advancing the fractional calculus in numerical modeling of engineering problems (Mehdinejadiani

*et al.*2013a, 2013b), in this study, the precision of the numerical simulation of pollution transmission in rivers is improved with the fractional approach. In other words, in this paper, the numerical model is developed based on the FRADE model. To validate the performance of the developed numerical model for theoretical and practical purposes, comparisons are conducted with the analytical solution of the MADE model and observed data from the Severn River in the UK.

## 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

*u*is the apparent flow velocity and

_{a}*D*is the corrected dispersion coefficient.

_{a}### Proposed numerical model for FRADE

*n*is the time and

*w*

_{j}

^{F}= (−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). In Equation (8), the first term on the right-hand side involves the Gaussian curve and

*FT*

_{1}is an attempt to model the fractional tail (

*FT*). In this study, the fully implicit scheme is developed and used as Equation (9): where

*m*= 1,2,3,. . . ,

*N*− 1 and

*j*=

*i*

*+*

*1-m*. Equation (9) can be rewritten as Equation (10): where and

*R*=

*C*

_{m}

^{n}+ .

*m*

*=*

*1*to

*N*−

*1*, 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.

. | Factor F. | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

Memory length j
. | 2.000 . | 1.900 . | 1.800 . | 1.700 . | 1.600 . | 1.500 . | 1.400 . | 1.300 . | 1.200 . | 1.100 . | 1.000 . |

0 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 |

1 | −2.000 | −1.90 | −1.80 | −1.70 | −1.60 | −1.50 | −1.40 | −1.30 | −1.20 | −1.10 | −1.00 |

2 | 1.000 | 0.855 | 0.720 | 0.595 | 0.480 | 0.375 | 0.280 | 0.195 | 0.120 | 0.055 | 0.000 |

3 | 0.000 | 0.029 | 0.048 | 0.060 | 0.064 | 0.063 | 0.056 | 0.046 | 0.032 | 0.017 | 0.000 |

4 | 0.000 | 0.008 | 0.014 | 0.019 | 0.022 | 0.023 | 0.022 | 0.019 | 0.014 | 0.008 | 0.000 |

5 | 0.000 | 0.003 | 0.006 | 0.009 | 0.011 | 0.012 | 0.012 | 0.010 | 0.008 | 0.005 | 0.000 |

6 | 0.000 | 0.002 | 0.003 | 0.005 | 0.006 | 0.007 | 0.007 | 0.006 | 0.005 | 0.003 | 0.000 |

7 | 0.000 | 0.001 | 0.002 | 0.003 | 0.004 | 0.004 | 0.005 | 0.004 | 0.004 | 0.002 | 0.000 |

8 | 0.000 | 0.001 | 0.001 | 0.002 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.002 | 0.000 |

9 | 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 j
. | 2.000 . | 1.900 . | 1.800 . | 1.700 . | 1.600 . | 1.500 . | 1.400 . | 1.300 . | 1.200 . | 1.100 . | 1.000 . |

0 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 |

1 | −2.000 | −1.90 | −1.80 | −1.70 | −1.60 | −1.50 | −1.40 | −1.30 | −1.20 | −1.10 | −1.00 |

2 | 1.000 | 0.855 | 0.720 | 0.595 | 0.480 | 0.375 | 0.280 | 0.195 | 0.120 | 0.055 | 0.000 |

3 | 0.000 | 0.029 | 0.048 | 0.060 | 0.064 | 0.063 | 0.056 | 0.046 | 0.032 | 0.017 | 0.000 |

4 | 0.000 | 0.008 | 0.014 | 0.019 | 0.022 | 0.023 | 0.022 | 0.019 | 0.014 | 0.008 | 0.000 |

5 | 0.000 | 0.003 | 0.006 | 0.009 | 0.011 | 0.012 | 0.012 | 0.010 | 0.008 | 0.005 | 0.000 |

6 | 0.000 | 0.002 | 0.003 | 0.005 | 0.006 | 0.007 | 0.007 | 0.006 | 0.005 | 0.003 | 0.000 |

7 | 0.000 | 0.001 | 0.002 | 0.003 | 0.004 | 0.004 | 0.005 | 0.004 | 0.004 | 0.002 | 0.000 |

8 | 0.000 | 0.001 | 0.001 | 0.002 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.002 | 0.000 |

9 | 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 |

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

*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.

## RESULTS AND DISCUSSION

^{2}/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.

*u*= 1.1 m/s and

*D*= 10.2 L

_{L}^{2}/s. As shown in this figure, by decreasing the concentration peak, the skewness of the curves increases.

*ɛ*= 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.

*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.

Station . | Spatial step (m) . | Time step (s) . | Fractional order (F factor)
. | Dispersion coefficient . |
---|---|---|---|---|

A | 4 | 2 | 1.68 | 59.2 |

B | 4 | 2 | 1.68 | 19.2 |

C | 4 | 2 | 1.68 | 12.8 |

D | 4 | 2 | 1.59 | 26.3 |

E | 4 | 2 | 1.53 | 37.4 |

F | 4 | 2 | 1.48 | 29.7 |

Station . | Spatial step (m) . | Time step (s) . | Fractional order (F factor)
. | Dispersion coefficient . |
---|---|---|---|---|

A | 4 | 2 | 1.68 | 59.2 |

B | 4 | 2 | 1.68 | 19.2 |

C | 4 | 2 | 1.68 | 12.8 |

D | 4 | 2 | 1.59 | 26.3 |

E | 4 | 2 | 1.53 | 37.4 |

F | 4 | 2 | 1.48 | 29.7 |

## 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.