Simulation of Diffusion Advection Phenomena with Retardation and Source/Sink Terms Using a Series of Completely Mixed Reactors

In this research, pollution transport in a 1-D domain with a source/sink and a retardation term is simulated using a series of completely mixed reactors. It is shown that as the volume of these reactors approaches zero, the governing equation of the system becomes identical to the advection diffusion partial differential equation. By application of mass balance to a finite number of completely mixed reactors, a set of ordinary differential equations which is equivalent to spatially discretized form of the equation is obtained in a completely natural and intuitive way. More importantly, the proposed innovative model can offer a conceptual insight into pollution transport phenomena and is flexible to any change for different problems.


Introduction
Pollution transport is one of the most significant problems in environmental engineering and mainly takes place as a result of advection, diffusion, and reactions in rivers. Simulation of such problems is mostly concerned with solving differential equations. Obtaining an analytical solution is not always possible, and therefore numerical methods are mostly used to solve these equations. On the other hand, a comprehensive model which is able to simulate the advection-diffusion process in a river physically, based on a familiar concept such as mass balance, may yield to establishment of a more welldefined scheme rather than those presented by numerical methods, typically. As physical sense can be highlighted in the model, implementing new initial and boundary conditions as well as the other desirable changes in the problem statement can be carried out in a more obvious framework.
Application of mass balance in completely mixed reactors for transport simulation has been used extensively by researchers (Hill and Root 2014;Carberry 2001;Jakobsen 2008;Mann 2009;Salmi et al. 2010;Davis and Davis 2012). Advection-diffusion simulation has been studied using feedforward and back-feed reactors by Chapra (1997). Dispersion effects have been addressed by Metcalf and Eddy (Burton et al. 2014). Effects of chemical reactions in pollutants fate using completely mixed reactors have been investigated by Hessel et al. (2004) and have been employed to indoor air quality simulations by Hayes (1991) and WHO (2015). As far as the authors are concerned, diffusion of pollution by 1 3 completely mixed reactors has not been simulated using flow rates to/from a reactor in feedforward and back-feed systems.
In this study, it is shown that application of mass balance to a series of reactors transforms the transport equation into a system of linear equations. This transformation is equivalent to application of numerical methods for solving the governing partial differential equation (PDE) of mass transport. Existence of source/sink and reactions has also been formulated. Furthermore, these formulations are verified using several examples and comparing the results with the analytical solution.

Model Construction
It is assumed that pollution transport in a 1-D river may be simulated by a series of completely mixed reactors. In this model, each reactor is connected to its left and right reactors via pipes. Water storage and reactions take place in these reactors, and driving forces for pollution transport (i.e., diffusion and advection) are provided by implementing pumps.

Diffusion Simulation
Diffusion is a type of random motion where molecular movements to left and right have the same probability. In order to model diffusion in the proposed model, an arbitrary reactor (e.g., reactor i, shown in Fig. 1a) is connected to its adjacent reactors at its left and right (reactors i − 1 and i + 1) by two identical pumps: i, i + 1. Using this concept, transmission of pollution to left and right with the same probability is allowed. This process is applied for all other reactors, as shown in Fig. 1a, to model diffusion process in the whole domain. As pumps start working simultaneously, pollution transports from time 0 onward. Fig. 1 a Simulation of diffusion using completely mixed reactors in series. b Simulation of advection using completely mixed reactors in series. c Simulation of diffusion-advection in a 1-D domain using completely mixed reactors

Advection Simulation
Advection is the transport of pollution by bulk flow in river. In order to model advection (form left to the right part of the system), each reactor is connected to its right reactor via a pipe. The driving force for advection is provided by implementing a pump at downstream (the most right reactor) having a flow rate of Q Adv (Fig. 1b). This layout allows the bulk transport of pollution from the left to the right of the system. Pollution concentration in the advective pipes is assumed to vary linearly.

Diffusion-Advection Simulation
By combination of the two mentioned systems, diffusion-advection may be modeled using a series of completely mixed reactors as shown in Fig. 1c.

Governing Equation
In order to derive the governing partial differential equation of the proposed model, it is assumed that the model consists of infinitely small reactors with a volume of dV. Since pollution varies linearly in the advective pipes, the concentration of pollution in these pipes is assumed to be the average of both their end concentrations. Mass balance equation for an interior reactor i shown in Fig. 2a is written as: (1) m in −ṁ out = dm dt .
Substituting QC for ṁ yields: Simplification of the above equations results in: (2) which shows that the governing equation for the proposed model is identical to the diffusion advection equation.

Discretization of the Governing Equation
In order to solve Eq. (7), one may use numerical methods to discretize this equation both spatially and temporally. In this research, however, a new procedure is followed for spatially discretizing Eq. (7). Writing mass balance equation for an arbitrary reactor (i), shown in Fig. 2b, yields: Simplification of Eq. (8) yields: Application of Eq. (9) to all interior nodes yields a system of linear first-order differential equations which leads to closure when considering the boundary nodes. The time derivative in Eq. (9) and similar equations may be evaluated using forward, backward or Crank-Nicolson finite difference schemes.

Boundary Nodes
In order to write mass balance equation for boundary reactors, two types of boundaries are investigated in this research: infinite source and impermeable boundary.
• Infinite source: when a boundary is an infinite source such as lakes, concentration in this boundary is specified and constant over time (Fig. 3a). In order to model this boundary, the reactor at that boundary is considered to have an infinite volume with a constant concentration. Inflow and outflow to and from this reactor do not change its volume nor its concentration. Application of mass balance to this reactor yields: As the volume of the boundary reactor approaches infinity, the left-hand side of the above equation becomes zero, reflecting the fact that the concentration of the reactor remains constant over time. • Impermeable boundary implies the existence of an impermeable barrier at the end of a river, such as a rock. Flow and pollution flux are zero across this boundary, and as a result there would be no advection in the domain. To model this boundary using complete mix reactors, the reactor at this boundary (reactor 1) is not connected to any exterior reactor, as shown in Fig. 3b, resulting in no outward flux conditions. Application of mass balance for the considered reactor results:

Forcing Function
This term is called source/sink in environmental modeling and may be a point source/sink or a distributed source/sink term. These terms imply any external pollution recharge or discharge to or from a reactor. In order to model a point source/sink, mass balance is simply written for that reactor (e.g., reactor i). For a distributed source/sink, the forcing function is integrated over the midpoint of the pipe i − 1, i

Reaction
This implies any mass generation or loss due to reactions depending on the concentration. The most common example is retardation pollution. Assuming a mass decay rate of kC n , mass balance equation for reactor i is written as: where k is the decay rate and n is an exponent indicating the order of reaction. (12)

Examples and Results
In order to illustrate how the proposed model may be employed for solving pollution transport and also to show its capabilities, several diffusion advection examples are solved.
3.1 In this example, at first advection diffusion in a domain with homogeneous Dirichlet boundary conditions and a point source at the center are considered: Simulation of this example using the proposed model is illustrated in Fig. 4.
In order to solve this example, the problem is discretized using 11 reactors (9 interior reactors and 2 boundary reactors). The reactors are assumed to have a unit volume, V = 1, and ∆x is considered to be 1. Application of Eqs. (9), (10), and (13) to interior reactors, boundary reactors, and the reactor with the source term, respectively, results: C 10 +C 11 2 = dC 11 dt Boundary reactors In Eq. (3), V 1 and V 11 are assumed to have an infinite volume. For numerical purposes, their volume is assumed to be 1000. Application of the above equations results in a firstorder system of linear differential equations. Using forward finite difference over time with a time step of ∆t = 0.01 s, the solutions may be obtained. The results for D = 1 m 2 /s and u = 0.0, 1.0, and 2.0 m/s are shown in Fig. 5.
It is observed that the results are in good agreement with the analytical solution. The plot shows that in the absence of advection the distribution of concentration is symmetric, indicating that there is a symmetric diffusion. The plot also shows that as time increases the concentration through the entire domain dampens. The plot also depicts that in case of advection, concentration is deviated to the right as the result of velocity.
As the 2nd part of this example, it is assumed that the pollution decays by a first-order reaction and a reaction constant of 0.1 s −1 , n = 1, k = 0.1 s −1 . Application of Eq. (15) for interior, boundary, and central reactors yields:  The results for D = 1 m 2 /s, u = 1 m/s and values of 0.1 and 0.5 s −1 for k are shown in Fig. 6 and show the capability of the model for considering reactions in diffusion advection problems. The root-mean-square error (RMSE) over the entire domain for this example (for the two parts) at different times is reported in Table 1. The results justify employing completely mixed reactors for advection diffusion simulation.
It is seen that when the retardation term increases, the concentration of pollution decreases throughout the entire domain.

3.2
The second problem considers an advection diffusion problem with a first-order decay reaction, a non-homogeneous boundary condition at the left end, and a homogeneous boundary condition at the other end. Since the pollution decays, the left boundary condition does not remain constant over time: In order to solve this problem, the same procedures described in the previous example are followed for all reactors except for the left one. This problem is shown in Fig. 7.
Since the concentration in the left boundary reactor decays linearly with time, mass balance application to this reactor results in: where V 1 is considered to be 1000 times Vi. The results of solving this problem using the proposed model are shown in Fig. 8 and given in Table 2.
The plot shows that as time increases from 0 to 3, concentration is pushed to the right side because of the velocity. It also indicates that as the retardation term increases, due to the degradation of the pollution, the concentration decreases.
As in the previous example, the results verify the application of the proposed model as a numerical method for modeling transport processes.

3.3
The third example considers diffusion problem in a domain where the right boundary condition is an impermeable boundary. Because of the no flux condition, no advection may be present. This problem is stated as: This simulation of this example is shown in Fig. 9.
The procedure for solving this problem is the same as the two previous problems except for the rightmost reactor which models the no flux boundary condition. As is shown in Fig. 3b, the reactor in the right side of the boundary is only connected to the reactor at its left, imitating the no flux boundary condition. Formulation of all reactors is the same as in the previous problems except for the rightmost reactor which is: The results of solving this problem are shown in Fig. 10 and given in Table 3.

Conclusion
In this research, it was shown how to simulate diffusion advection transport with a retardation and a source/sink term using a series of completely mixed reactors in a 1-D domain such as rivers. By application of mass balance for each reactor, a system of first-order linear differential equations with respect to time was obtained. The system which was a spatially discretized form of the PDE governing pollution transport may then be solved to obtain the results of such simulations. In order to show the capabilities of the proposed model, three different examples were solved and root-mean-square errors were obtained. The results of the three examples showed that the solution obtained by the proposed model is in a good agreement with the analytical solutions. It can be confirmed that pollution transport in a 1-D domain including a retardation and/or a sink/source can be modeled using the proposed scheme accurately, along with a high degree of flexibility in imposing different conditions with more aspects of physical sense.

Appendix 1
In order to solve the three given examples in this paper, it should be noted that the diffusion advection equation with a linear retardation term can be transformed into the standard form of diffusion equation as follows: Solution for the first example: the original equation is transformed to: The eigenfunction for this problem is X(x) = sin n 10 x which is used for the expansion of other functions.

3
And w(x, t) is expanded as: And the solution is found as: By considering special times and points, these analytic results are given: w(x, t) = ∞ ∑ n=1 a n (t) sin n 10 x →̇a n (t) + D n 10 2 a n (t) x .

Appendix 2
Solution for the second example: Using the given transformation, this equation is transformed to: The eigenfunction for this problem is X(x) = sin n 10 x which is used for the expansion of other functions.
Since the problem has a non-homogeneous boundary condition, the following transformation may be used: By assuming a function that satisfies the non-homogeneous boundary conditions for s(x, t) such as: And substituting r(x, t) + s(x, t) for w(x, t), the following equation is achieved: Since s(x, t) is known, the two terms in the right-hand side of the above equation are expanded as: By considering special times and points, these analytic results are given: w(x, t) = r(x, t) + s(x, t).