Abstract
Two mathematical models are used to simulate water quality in a nonuniform flow stream. The first model is the hydrodynamic model that provides the velocity field and the elevation of water. The second model is the dispersion model that provides the pollutant concentration field. Both models are formulated in one-dimensional equations. The traditional Crank-Nicolson method is also used in the hydrodynamic model. At each step, the flow velocity fields calculated from the first model are the input into the second model as the field data. A modified MacCormack method is subsequently employed in the second model. This paper proposes a simply remarkable alteration to the MacCormack method so as to make it more accurate without any significant loss of computational efficiency. The results obtained indicate that the proposed modified MacCormack scheme does improve the prediction accuracy compared to that of the traditional MacCormack method.
1. Introduction
In general, the amount of pollution levels in a stream can be measured via data collection from a real of field data site. It is somehow rather difficult and complex, and the results obtained tentatively deviate in the measurement from one point in each time/place to another when the water flow in the stream is not uniform. In water quality modelling for nonuniform flow stream, the general governing equations used are the hydrodynamic model and the dispersion model. The one-dimensional shallow water equation and advection-dispersion-reaction equation is govern the first and the second models, respectively.
Numerous numerical techniques for solving such models are available. In [1], 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-dispersion-reaction equation, are presented in [2–6].
The nonuniform flow model requires data concerned with the 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 [7–10], they used the hydrodynamics model and the advection-dispersion equation to approximate the velocity of the water current in bay, uniform reservoir, and stream, respectively. Among these numerical techniques, the finite difference methods, including both explicit and implicit schemes, are mostly used for one-dimensional domain such as in longitudinal stream systems [11, 12].
There are two mathematical models used to simulate water quality in a nonuniform water flow systems. The first is the hydrodynamic model that provides the velocity field and the elevation of water. The second is the dispersion model that gives the pollutant concentration field. A couple of models are formulated in one-dimensional equations. The traditional Crank-Nicolson method is used in the hydrodynamic model. At each step, the calculated flow velocity fields of the first model are input into the second model as the field data [9, 10, 13].
The numerical techniques to solve the nonuniform flow of stream water quality model containing one-dimensional advection-dispersion-reaction equation have been presented in [10] using the fully implicit scheme: Crank-Nicolson method is used to solve the hydrodynamic model and backward time central space (BTCS) for dispersion model, respectively. In [13], the Crank-Nicolson method is also used to solve the hydrodynamic model, while the explicit Saulyev scheme is used to solve the dispersion model.
Their research on finite difference techniques for the dispersion model has concentrated on computation accuracy and numerical stability. Many complicated numerical techniques, such as the QUICK scheme, the Lax-Wendroff scheme, and the Crandall scheme, have been studied to increase performances. These techniques have focused on advantages in terms of stability and higher order accuracy [3].
The simple finite difference schemes become more attractive for model use. The simple explicit methods include the forward time-central space (FTCS) scheme, the MacCormack scheme, and the Saulyev scheme, and the implicit schemes include the backward time-central space (BTCS) scheme, and the Crank-Nicolson scheme [12]. These schemes are either first-order or second-order accurate and have the advantages in programming and computing without losing much accuracy and thus they are used for many model applications [3].
A third-order upwind scheme for the advection-diffusion equation using a simple spreadsheets simulation is proposed in [14]. In [15], a new flux splitting scheme is proposed. The scheme is robust and converges as fast as the Roe Splitting. The Godunov-mixed methods for advection-dispersion equations are introduced in [16]. A time-splitting approach for advection-dispersion equations is also considered. In addition, [17] proposes a time-splitting method for multidimensional advection-diffusion equations that advection is approximated by a Godunov-type procedure, and diffusion is approximated by a low-order mixed finite element method. In [18], the flux-limiting solution techniques for simulation of reaction diffusion convection system are proposed. A composite scheme to solves the scalar transport equation in a two-dimensional space that accurately resolve sharp profiles in the flow is introduced. The total variation diminishing implicit Runge-Kutta methods for dissipative advection-diffusion problems in astrophysics is proposed in [19]. They derive dissipative space discretizations and demonstrate that together with specially adapted total-variation-diminishing (TVD) or strongly stable Runge-Kutta time discretizations with adaptive step-size control this yields reliable and efficient integrators for the underlying high-dimensional nonlinear evolution equations.
In this research, we propose simple revisions to the MacCormack scheme that improve its accuracy for the problem of water quality measurement in a nonuniform water flow in a stream. In the following sections, the formulation of the traditional MacCormack scheme is reviewed. The revision of the modified MacCormack scheme is proposed.
The results from the hydrodynamic model are the data of the water flow velocity for the advection-dispersion-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 guarantees the accuracy of the approximate solution is presented in [9, 10, 13].
The stream has a simple one-space dimension as shown in Figure 1. By averaging the equation over the depth, discarding the term due to the Coriolis force, it follows that the one-dimensional shallow water and the advection-dispersion-reaction equations are applicable. We use the Crank-Nicolson scheme, the traditional MacCormack scheme, and the Modified MacCormack scheme to approximate the velocity, the elevation, and the pollutant concentration from the first and the second models, respectively.
2. Model Formulation
2.1. The Hydrodynamic Model
In this section, we derive a simple hydrodynamic model for describing water current and elevation by one-dimensional shallow water equation. We make the usual assumption in the continuity and momentum balance; that is, we assume that the Coriolis, shearing stresses, and the surface wind are small [7, 9, 10, 20]; we obtain the one-dimensional shallow water equations: where is the longitudinal distance along the stream (m), is time (s), is the depth measured from the mean water level to the stream bed (m), is the elevation from the mean water level to the temporary water surface or the tidal elevation (m/s), and is the velocity components (m/s), for all .
Assume that is a constant and . Then (1) lead to
We will consider the equation in the dimensionless problem by letting , and . Substituting them into (2) leads to
In [9, 10, 13], they introduce a damping term into (3) to represent the frictional forces due to the drag of sides of the stream, with the initial conditions at and are and . The boundary conditions for are specified: at and at . The system of (4) is called the damped equation. We solve the damped equation by using the finite difference method. In order to solve (4) in , it is convenient to use for and , respectively: with the initial conditions , at , and the boundary conditions and at .
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 [2–4, 6, 10] as follows: where is the concentration averaged in depth at the point and at time , is the diffusion coefficient, is the mass decay rate, and is the velocity component for all . We will consider the model with the following conditions. The initial condition at for all . The boundary conditions at and at where is a constant.
3. Crank-Nicolson Method for the Hydrodynamic Model
The hydrodynamic model provides the velocity field and elevation of the water. Then the calculated results of the model will be the input into the dispersion model which provides the pollutant concentration field. We will follow the numerical techniques of [9]. To find the water velocity and water elevation from (5), we make the following change to variables and substitute it into (5). Therefore, Equations (7) can be written in the matrix form That is where with the initial condition at . The left boundary condition for is specified: and , and the right boundary condition for is specified: and .
We now discretize (9) by dividing the interval into subintervals such that and the interval into subintervals such that . We can then approximate by , value of the difference approximation of at point and , where and , and similarly defined for and . The grid point is defined by for all and for all in which and are positive integers. Using the Crank-Nicolson method [21] to (9), the following finite difference equation can be obtained: where is the unit matrix of order 2, and . Applying the initial and boundary conditions given in (7), it can be obtained the general form where
where , and for all . The Crank-Nicolson scheme is unconditionally stable [12, 21].
4. A Modified MacCormack Scheme for the Advection-Dispersion-Reaction Equation
4.1. The Traditional MacCormack Scheme
First of all, we consider the traditional MacCormack scheme. The scheme is an explicit finite difference scheme with predictor-corrector two-step method. The first step is a modification of forward time central space (FTCS) by changing the central space evaluation at time to a forward space evaluation. This step is a forward time forward space (FTFS) scheme. The FTFS scheme approximates the temporal and spacial derivatives and the decay in (6) with the following discretization.
We can then approximate by , the value of the difference approximation of at point and , where and . The grid point is defined by for all and for all in which and are positive integers. Taking the forward time forward space technique [3, 21] into (6), we get the following discretization:
Note that are obtained by the Crank-Nicolson method with the hydrodynamic model of (5) that are presented in [9, 10, 13].
Substituting (15) into (6), we get for and . Substitute the difference equation into (16), and then define slope as Let and , and then define and . Equation (17) takes a simplified form: or
For upper boundary, where , plug the known value of the left boundary to (19) in the right hand side; we obtain
For the lower boundary, where , substitute the approximate unknown value of the right boundary by the forward difference approximation to . Let and rearrange; we obtain
Taking the Euler formula, we obtain the MacCormack predictor step formulation:
The second step is a modified backward time central space (BTCS) scheme by changing the central space evaluation time with a backward space evaluation. It is essentially a backward time backward space (BTBS) scheme. The BTBS scheme approximates the temporal and spacial derivatives and the decay in (6) with the following discretization:
Because the values at time level have been calculated in predictor step, the second step is also explicit. It follows that the slope base on their predictor points can be calculated as follows:
For upper boundary, where , plug the known value of the left boundary to (24) in the right hand side. We obtain
For the lower boundary, where , substitute the approximate unknown value of the right boundary by the backward difference approximation to . Let and rearrange; then, we obtain
From both of the steps, the MacCormack scheme takes the following form:
The MacCormack scheme is conditionally stable subject to constraints in (16). The stability requirements for the scheme are [22] where is the diffusion number (dimensionless) and is the advection number or Courant number (dimensionless).
4.2. The Modified MacCormack Scheme
Since the derivative approximation during discretization is not centered, numerical dispersion will be introduced. The dispersion coefficients used in the dispersion model would take the value obtained by subtracting the numerical dispersion from the real data of the stream. The amounts of the numerical dispersion introduced by backward space denoted by and forward time denoted by schemes as follow [3, 12]:
There are temporal and spacial numerical dispersion in both predictor and corrector steps since the scheme uses forward time forward space difference for prediction and backward time backward space difference for correction. From (29), the numerical dispersion for forward time forward space prediction step and backward time backward space correction step are
The modified MacCormack scheme uses the following corrected dispersion, rather than the real dispersion coefficients for calculation in both prediction and correction steps: where is the dispersion coefficient used in the prediction step and is the dispersion coefficient used in the correction step.
The modified MacCormack scheme is conditionally stable subject to the constraint in (16). The stability requirements for the scheme: where the maximum of numerical dispersion coefficients is .
5. The Accuracy of the Hydrodynamic Approximation
It is not hard to find the analytical solution in (5) with . By changing of variables, and for some by substituting it in (5). Using a separable variables technique, we can obtain a solution [10]: where and . However, it is not easy to find the analytical solution of (5). We use the solution obtained in (33) to verify our approximate solution obtained by the Crank-Nicolson method equation (13). Actually when the Crank-Nicolson method is used, we get the approximate solution both and . We assume that when we get a good approximation for , this implies that the method gives a good approximation for . The verification of the approximate solution is shown in Figure 2.
Figure 2 shows the comparison between the analytical solutions and the approximate solutions only at the end of the domain.
Unfortunately, the analytical solutions of the hydrodynamic model could not be found over the entire domain [10]. This implied that the analytical solutions of the dispersion model could not be computed at any points in the domain as well.
6. Application to the Stream Water Quality Assessment Problem
Suppose that the measurement of pollutant concentration in a nonuniform flow stream 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 the discharge point is (mg/L) at for all and (mg/L) at . The elevation of water at the discharge point can be described as a function (m) for all , and the elevation does not change at (km) The physical parameters of the stream system are diffusion coefficient (m2/s) and a first-order reaction rate s−1. In the analysis conducted in this study, meshing the stream into 40 elements with , and the time increment is 0.4 (s) with , characterizing a one-dimensional flow. Using the Crank-Nicolson method of [9, 10, 13], it can be obtained the water velocity in Table 1 and Figure 3. Next, the approximate water velocity can be plugged into the traditional MacCormack scheme in (27). We also plug the approximate water velocity into the modified MacCormack scheme (27) with numerical dispersion coefficients (31). The approximation of pollutant concentrations of both schemes is shown in Tables 2 and 3 and Figures 4 and 5. The comparison of traditional MacCormack and modified MacCormack is shown in Figure 6.
(a)
(b)
(c)
(d)
7. Discussion and Conclusions
The approximation of the pollutant concentrations of the traditional and modified MacCormack schemes is shown in Tables 2 and 3. The real-world problems require a small amount of time interval in obtaining accurate solutions. Unfortunately, the analytical solutions of the hydrodynamic model could not be found over the entire domain. This also implies that the analytical solutions of dispersion model could not work out at any point on the entire domain as well [13].
In [13], it is revealed that the diffusion coefficients of the pollutant matter can reduce the concentration of the 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.
We propose a modified MacCormack scheme by adding a simple revision to the traditional MacCormack scheme. The numerical dispersion has been introduced because the derivative approximation during discretization is not centered. The traditional MacCormack scheme shows excessive dispersion effects for large time and space step lengths, significantly decreasing the efficiency of the traditional MacCormack scheme [3]. To eliminate the numerical dispersion effect, the modified MacCormack scheme for the nonuniform flow is proposed. Though revision shows a good agreement in accuracy with the original one, the modified MacCormack scheme becomes less efficient than the traditional MacCormack scheme.
In this paper, the hydrodynamic model and the convection-diffusion-reaction equation can be combined to approximate the pollutant concentration in a stream when the current reflecting 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 discharge point, can be obtained. Both of the traditional and the modified MacCormack schemes can be used in the dispersion model since the scheme is very simple to implement. By both of the traditional and the modified MacCormack finite difference formulations, we obtain that the proposed technique is applicable and economical to be used in the real-world problem due to its simplicity to program and the straight forwardness of the implementation. It is also possible to find tentative better locations and better periods of time of the different discharge points to a stream.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The author greatly appreciates all valuable comments received from Prof. Chatchai Leenawong and the referees.