- Email: [email protected]
Similarity of kinematic and diffusive waves: a comment on accuracy criteria for linearised diffusion wave flood routing
Journal of Hydrology 216 (1999) 248–251
Similarity of kinematic and diffusive waves: a comment on accuracy criteria for linearised diffusion wave flood routing. By K. Bajracharya and D.A. Barry [Journal of Hydrology, vol. 195 (1997), 200-217] R. Szymkiewicz Department of Hydraulics, Technical University of Gdan´sk 80-952Gdan´sk, Poland
According to the title, in the referred paper an accuracy analysis for linearised diffusion flood routing should be presented. The main purpose of this analysis was to obtain the proper values of mesh dimensions as Dt and Dx ensuring acceptable agreement between the results of computation and observed data. However, the Authors presented an accuracy analysis for kinematic wave solved by a box scheme containing two weighting parameters. According to them the presented analysis concerns diffusive wave equation. The Authors accepted the solution of kinematic wave equation containing large numerical error as a solution of diffusive wave equation. In fact the results of calculations can be similar but it does not mean that the equations giving them are also similar. For this reason I cannot agree with the title of referred paper. It contains an analysis concerning kinematic wave equation or more generally an advection equation, which has nothing in common with the accuracy analysis of numerical solution of advection-diffusion equation as diffusive wave model. Consequently the presented conclusions concern the kinematic wave equation. Additionally, some of them are well known and can be formulated without any numerical calculations because they result from theory of numerical methods for partial differential equations. The following are the reasons which allow me to formulate this opinion.
First let me recollect some information about the Muskingum model, which is a lumped model in the form of a linear ordinary differential equation containing two parameters: X and K. The parameter K represents the wave travel time between two river sections whereas the parameters X (called by Authors u ) has no physical interpretation. Cunge (1969) showed a kinematic character of Muskingum model. Moreover he stated that the wave attenuation in Muskingum model is caused by numerical diffusion. When the trapezodial rule (v 1/2) is used to time integration, the numerical diffusion is controlled by u and D x only. Cunge proposed to accept the value u producing numerical diffusion equivalent to the hydraulic one existing in diffusive wave model. This form of Muskingum model is called Muskingum-Cunge model. It is possible to show readily that the Muskingum model is a semidiscrete form of linear kinematic wave equation obtained by space derivative approximation at point between nodes xi⫺1 and xi. Its location is defined by parameter u . The approximation is of II order accuracy with regard to x for u 1/2 only. For other values of u it is an approximation of I order accuracy. As the result of poor spatial approximation dissipation error in the form of numerical diffusion occurs. It ensures a flood wave attenuation by kinematic wave model. One can add that u ⱕ 1/2 is
0022-1694/99/$ - see front matter 䉷 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(98)00276-5
R. Szymkiewicz / Journal of Hydrology 216 (1999) 248–251
necessary to ensure numerical stability of solution. Cunge (1967) has also presented these properties. The Authors use the term ‘‘accuracy’’ in two meanings. The first meaning corresponds to the one used in mathematics when the numerical solution is related with exact solution. The second meaning deals with the relation between numerical solution and observations. I have to add that the Authors also use the term ‘‘optimal solution’’ which, in my opinion, is apt. Those terms have different sense and they should be distinguished. The solution similar to the observations can be obtained even when the applied model represents poor accuracy of approximation. This statement seems to be a paradox but in hydrology it is sometimes true. An example of this is the Muskingum model as well as the linear and non-linear reservoir models. Good agreement between the calculations and observations is ensured by numerical error produced by applied scheme. We have to remember that in reality an accuracy analysis of Muskingum method is an accuracy of kinematic wave because it is governing equation for Muskingum model. However, the Authors first on page 200 say that the MuskingumCunge model is a particular case of difference approximation of kinematic wave and then on page 201 state that this same model is now an approximation of diffusive wave equation. It seems that the Muskingum model cannot be an approximation of kinematic and diffusive waves simultaneously. They are completely different differential equations having different exact solutions and additionally they are solved by different numerical schemes. Similarity of exact solution of diffusive wave with numerical solution of kinematic wave ensured by large numerical error does not allow to formulate the conclusion about a similarity of equations. This problem is explained by convergence condition, which should be satisfied by each numerical scheme. As was stated before, it is easy to show that the Muskingum method is an approximation of kinematic wave, even from Eq. 7, which is a modified governing equation, results that it is consistent with kinematic wave but not with diffusive wave. For Dx ! 0 and Dt ! 0 it does not tend to Eq. 1 but tends to Eq. 2. Because simultaneously a numerical solution is stable then according to the Lax theorem the convergence is ensured. However, the Muskingum solution converges to solution of linear kinematic wave what
249
can be showed by numerical experiments for decreasing mesh dimensions. The final conclusion is as follows: Eq. (7) cannot serve the accuracy analysis of diffusive wave. Because Eq. (5a) approximates a kinematic wave equation, the Authors’ analysis presented in the paper concerns this type of equation which of course coincides with analysis of Muskingum model. The scheme (5a) is an approximation of I order accuracy of kinematic wave as long as the coefficient of numerical diffusion Dn is different to zero. It is defined as follows: Dn
CDx ÿ 1 ⫺ 2u ⫹ 2v ⫺ 1Cr 2
A
where C is the kinematic velocity; Dx the distance step; u ,v are weighting parameters; Cr the advective Courant number.Therefore when Dn 苷 0 one cannot say about an accuracy of II, III or IV order as it is presented by Authors in chapters 2.3, 2.4, 2.5. It is still an accuracy of I order even if the terms containing higher order derivatives disappear. In this case the cancelling of terms with III or IV order derivatives has nothing common with approximation accuracy of III or IV order not only for diffusive wave but as well for kinematic wave. One can add that cancelling of higher terms does not change the results considerably when the term with II order derivative is saved. This fact can be easily confirmed by solving Eq. (5a). With v 0 and u 0 one obtains the following relation: CDx j j j Q Q ⫺ ⫺ Q B Qj⫹1 i j i i⫺1 Dx which coincides with Eq. (15): ÿ 1 ⫺ Cr Qji ⫹ Cr Qji⫺1 Qj⫹1 i
C
It is a very well known and examined upwind scheme, widely applied in fluid mechanics to present the problem of numerical diffusion. It is described in may books for instance by Fletcher (1991),Abbott and Basco (1989). For this reason, I do not know why Authors call it ‘‘ mixing cell model’’ additionally derived in a very complicated manner. One can add that the scheme (C) coincides with linear reservoir model with the explicit Euler method of integration in time. It is known that even if the discontinuities exist in initial or boundary conditions this scheme
250
R. Szymkiewicz / Journal of Hydrology 216 (1999) 248–251
always ensures solution without any oscillations but strongly smoothed by numerical diffusion. Additionally the peak wave is always in-phase with the true solution (Abbott and Basco, 1989). The numerical diffusion increases when Courant number decreases. For Cr 1 an exact solution of kinematic wave is obtained. The exact solution is given by scheme when in modified Eq. (7) all terms on the right side disappear. However, for v 0, u 0 and Cr 1/2 the term with III order derivative is cancelled but the obtained solution is not more accurate. It is evident because the truncation error is usually dominated by its first term, which in this case is diffusive one. If the main purpose of referred paper was to obtain an optimal solution of linear diffusive wave using general difference implicit scheme, the solution of this problem seems relatively simple. One can use Crank-Nicolson scheme applied by authors to reproduce ‘‘an exact solution’’. It is well known that Crank-Nicolson scheme being an approximation of II order with regard to x as well to t is numerical dissipation-free. It does not produce any numerical diffusion but unfortunately it always generates dispersion error. This causes that Crank-Nicolson scheme is useless for pure advection equation. Whereas it can be very useful for advection-diffusion equation when the Peclet number is not great than 2: P
CDx ⱕ2 D
D
This condition ensures a smooth solution. Therefore when we accept: Dx ⱕ
2D C
E
a satisfying accuracy of solution can be expected. For a diffusive wave moving in rectangular channel one obtains: C mU D
Q UH 2Bs 2s
F G
Where m 5/3 (for Manning formula); U–average flow velocity; H–depth; B–channel width; s–bed slope.With these relations inequality (E) takes from Dx ⱕ
H ms
H
For real data this condition can be readily satisfied and Crank-Nicolson scheme will ensure success practically always. The numerical solution of diffusive wave by this scheme will tend to exact one when Dx ! 0 and Dt ! 0 because Crank-Nicolson scheme satisfies the consistence condition and is absolutely stable. Moreover relatively small computational effort is needed to obtain numerical solution. Nevertheless, the Authors as well as many hydrologists prefer to reach similar results by using the Muskingum model. By doing this they face some difficulties as for obtaining a good agreement between the calculations and observations, and to attain this a proper choice of Dx and Dtis necessary. It should be remembered that the purpose of this choice is to produce sufficiently large dissipation error, determined by Dx and Dt which ensures an attenuation of flood wave as diffusive wave model. For this reason the readers who are familiar with the theory of numerical methods may be astonished by the Author’s statement, that the increase of solution accuracy is caused by the increase in Dx Moreover the excellent results ensure the less accurate scheme (v 0, u 0). This is, without any doubt, the serious weakness of the hydrological flood routing models such as the Muskingum model. In these models the effect of lacking physical dissipation is reproduced by numerical error. This fact should encourage the hydrologists to reflect about the usefulness of application of this kind of models. On the contrary, diffusive wave model contains physically founded hydraulic diffusion mechanism and the numerical methods applied to solve it are correct. When we use the Muskingum model in the form of Eq. 5a, the accepted values of v and u should ensure a relatively small value of Dx defined by Eq. 8. In an other case, when u ! 1=2 and v ! 1=2 the approximation accuracy of kinematic wave increases. According to Eq. (8) we observe that to obtain a numerical diffusion comparable with hydraulic one in diffusive wave, Dx ! ∞. However, from the theory of numerical methods it is known that the numerical grid allows to reproduce the waves longer than 2Dx only. Generally speaking, the different methods lead to a solution that represents the long wave properties only, whereas the short waves are considerably disturbed. Therefore if the flood wave has large gradients, the oscillations of the solution should be
R. Szymkiewicz / Journal of Hydrology 216 (1999) 248–251
expected even when Dx satisfies Eq. (8). For this reason it is rather impossible to obtain some general formulas defining the optimal values of distance step Dx and time step Dt for Muskingum model because they can depend on the character of flood wave also. These problems are connected with a general principle; the increasing order of approximation accuracy does not have to improve the accuracy of solution when the grid is sparse. Inversely, the solution accuracy suffers. The remarkable improvement of accuracy of the solution occurs when the increasing order of approximation appears together with decreasing mesh dimensions (Fletcher, 1991). Unfortunately it is not interesting for users of hydrological models because this process, reducing the numerical error
251
leads to the exact solution of kinematic wave which, as is known, for C const corresponds to pure translation of flood wave without any shape deformation.
References Abbott, M.B., Basco, D.R., 1989. Computational Fluid Dynamics. Longman Scientific and Technical. Cunge, J., 1969. On the subject of a flood propagation computation method (Muskingum method). J. of Hydraulic Research 7, 205– 230. Fletcher, C.A., 1991. Computational Techniques for Fluid Dynamics, 1, Springer-Verlag, Germany.
Similarity of kinematic and diffusive waves: a comment on accuracy criteria for linearised diffusion wave flood routing. By K. Bajracharya and D.A. Barry [Journal of Hydrology, vol. 195 (1997), 200-217] R. Szymkiewicz Department of Hydraulics, Technical University of Gdan´sk 80-952Gdan´sk, Poland
According to the title, in the referred paper an accuracy analysis for linearised diffusion flood routing should be presented. The main purpose of this analysis was to obtain the proper values of mesh dimensions as Dt and Dx ensuring acceptable agreement between the results of computation and observed data. However, the Authors presented an accuracy analysis for kinematic wave solved by a box scheme containing two weighting parameters. According to them the presented analysis concerns diffusive wave equation. The Authors accepted the solution of kinematic wave equation containing large numerical error as a solution of diffusive wave equation. In fact the results of calculations can be similar but it does not mean that the equations giving them are also similar. For this reason I cannot agree with the title of referred paper. It contains an analysis concerning kinematic wave equation or more generally an advection equation, which has nothing in common with the accuracy analysis of numerical solution of advection-diffusion equation as diffusive wave model. Consequently the presented conclusions concern the kinematic wave equation. Additionally, some of them are well known and can be formulated without any numerical calculations because they result from theory of numerical methods for partial differential equations. The following are the reasons which allow me to formulate this opinion.
First let me recollect some information about the Muskingum model, which is a lumped model in the form of a linear ordinary differential equation containing two parameters: X and K. The parameter K represents the wave travel time between two river sections whereas the parameters X (called by Authors u ) has no physical interpretation. Cunge (1969) showed a kinematic character of Muskingum model. Moreover he stated that the wave attenuation in Muskingum model is caused by numerical diffusion. When the trapezodial rule (v 1/2) is used to time integration, the numerical diffusion is controlled by u and D x only. Cunge proposed to accept the value u producing numerical diffusion equivalent to the hydraulic one existing in diffusive wave model. This form of Muskingum model is called Muskingum-Cunge model. It is possible to show readily that the Muskingum model is a semidiscrete form of linear kinematic wave equation obtained by space derivative approximation at point between nodes xi⫺1 and xi. Its location is defined by parameter u . The approximation is of II order accuracy with regard to x for u 1/2 only. For other values of u it is an approximation of I order accuracy. As the result of poor spatial approximation dissipation error in the form of numerical diffusion occurs. It ensures a flood wave attenuation by kinematic wave model. One can add that u ⱕ 1/2 is
0022-1694/99/$ - see front matter 䉷 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(98)00276-5
R. Szymkiewicz / Journal of Hydrology 216 (1999) 248–251
necessary to ensure numerical stability of solution. Cunge (1967) has also presented these properties. The Authors use the term ‘‘accuracy’’ in two meanings. The first meaning corresponds to the one used in mathematics when the numerical solution is related with exact solution. The second meaning deals with the relation between numerical solution and observations. I have to add that the Authors also use the term ‘‘optimal solution’’ which, in my opinion, is apt. Those terms have different sense and they should be distinguished. The solution similar to the observations can be obtained even when the applied model represents poor accuracy of approximation. This statement seems to be a paradox but in hydrology it is sometimes true. An example of this is the Muskingum model as well as the linear and non-linear reservoir models. Good agreement between the calculations and observations is ensured by numerical error produced by applied scheme. We have to remember that in reality an accuracy analysis of Muskingum method is an accuracy of kinematic wave because it is governing equation for Muskingum model. However, the Authors first on page 200 say that the MuskingumCunge model is a particular case of difference approximation of kinematic wave and then on page 201 state that this same model is now an approximation of diffusive wave equation. It seems that the Muskingum model cannot be an approximation of kinematic and diffusive waves simultaneously. They are completely different differential equations having different exact solutions and additionally they are solved by different numerical schemes. Similarity of exact solution of diffusive wave with numerical solution of kinematic wave ensured by large numerical error does not allow to formulate the conclusion about a similarity of equations. This problem is explained by convergence condition, which should be satisfied by each numerical scheme. As was stated before, it is easy to show that the Muskingum method is an approximation of kinematic wave, even from Eq. 7, which is a modified governing equation, results that it is consistent with kinematic wave but not with diffusive wave. For Dx ! 0 and Dt ! 0 it does not tend to Eq. 1 but tends to Eq. 2. Because simultaneously a numerical solution is stable then according to the Lax theorem the convergence is ensured. However, the Muskingum solution converges to solution of linear kinematic wave what
249
can be showed by numerical experiments for decreasing mesh dimensions. The final conclusion is as follows: Eq. (7) cannot serve the accuracy analysis of diffusive wave. Because Eq. (5a) approximates a kinematic wave equation, the Authors’ analysis presented in the paper concerns this type of equation which of course coincides with analysis of Muskingum model. The scheme (5a) is an approximation of I order accuracy of kinematic wave as long as the coefficient of numerical diffusion Dn is different to zero. It is defined as follows: Dn
CDx ÿ 1 ⫺ 2u ⫹ 2v ⫺ 1Cr 2
A
where C is the kinematic velocity; Dx the distance step; u ,v are weighting parameters; Cr the advective Courant number.Therefore when Dn 苷 0 one cannot say about an accuracy of II, III or IV order as it is presented by Authors in chapters 2.3, 2.4, 2.5. It is still an accuracy of I order even if the terms containing higher order derivatives disappear. In this case the cancelling of terms with III or IV order derivatives has nothing common with approximation accuracy of III or IV order not only for diffusive wave but as well for kinematic wave. One can add that cancelling of higher terms does not change the results considerably when the term with II order derivative is saved. This fact can be easily confirmed by solving Eq. (5a). With v 0 and u 0 one obtains the following relation: CDx j j j Q Q ⫺ ⫺ Q B Qj⫹1 i j i i⫺1 Dx which coincides with Eq. (15): ÿ 1 ⫺ Cr Qji ⫹ Cr Qji⫺1 Qj⫹1 i
C
It is a very well known and examined upwind scheme, widely applied in fluid mechanics to present the problem of numerical diffusion. It is described in may books for instance by Fletcher (1991),Abbott and Basco (1989). For this reason, I do not know why Authors call it ‘‘ mixing cell model’’ additionally derived in a very complicated manner. One can add that the scheme (C) coincides with linear reservoir model with the explicit Euler method of integration in time. It is known that even if the discontinuities exist in initial or boundary conditions this scheme
250
R. Szymkiewicz / Journal of Hydrology 216 (1999) 248–251
always ensures solution without any oscillations but strongly smoothed by numerical diffusion. Additionally the peak wave is always in-phase with the true solution (Abbott and Basco, 1989). The numerical diffusion increases when Courant number decreases. For Cr 1 an exact solution of kinematic wave is obtained. The exact solution is given by scheme when in modified Eq. (7) all terms on the right side disappear. However, for v 0, u 0 and Cr 1/2 the term with III order derivative is cancelled but the obtained solution is not more accurate. It is evident because the truncation error is usually dominated by its first term, which in this case is diffusive one. If the main purpose of referred paper was to obtain an optimal solution of linear diffusive wave using general difference implicit scheme, the solution of this problem seems relatively simple. One can use Crank-Nicolson scheme applied by authors to reproduce ‘‘an exact solution’’. It is well known that Crank-Nicolson scheme being an approximation of II order with regard to x as well to t is numerical dissipation-free. It does not produce any numerical diffusion but unfortunately it always generates dispersion error. This causes that Crank-Nicolson scheme is useless for pure advection equation. Whereas it can be very useful for advection-diffusion equation when the Peclet number is not great than 2: P
CDx ⱕ2 D
D
This condition ensures a smooth solution. Therefore when we accept: Dx ⱕ
2D C
E
a satisfying accuracy of solution can be expected. For a diffusive wave moving in rectangular channel one obtains: C mU D
Q UH 2Bs 2s
F G
Where m 5/3 (for Manning formula); U–average flow velocity; H–depth; B–channel width; s–bed slope.With these relations inequality (E) takes from Dx ⱕ
H ms
H
For real data this condition can be readily satisfied and Crank-Nicolson scheme will ensure success practically always. The numerical solution of diffusive wave by this scheme will tend to exact one when Dx ! 0 and Dt ! 0 because Crank-Nicolson scheme satisfies the consistence condition and is absolutely stable. Moreover relatively small computational effort is needed to obtain numerical solution. Nevertheless, the Authors as well as many hydrologists prefer to reach similar results by using the Muskingum model. By doing this they face some difficulties as for obtaining a good agreement between the calculations and observations, and to attain this a proper choice of Dx and Dtis necessary. It should be remembered that the purpose of this choice is to produce sufficiently large dissipation error, determined by Dx and Dt which ensures an attenuation of flood wave as diffusive wave model. For this reason the readers who are familiar with the theory of numerical methods may be astonished by the Author’s statement, that the increase of solution accuracy is caused by the increase in Dx Moreover the excellent results ensure the less accurate scheme (v 0, u 0). This is, without any doubt, the serious weakness of the hydrological flood routing models such as the Muskingum model. In these models the effect of lacking physical dissipation is reproduced by numerical error. This fact should encourage the hydrologists to reflect about the usefulness of application of this kind of models. On the contrary, diffusive wave model contains physically founded hydraulic diffusion mechanism and the numerical methods applied to solve it are correct. When we use the Muskingum model in the form of Eq. 5a, the accepted values of v and u should ensure a relatively small value of Dx defined by Eq. 8. In an other case, when u ! 1=2 and v ! 1=2 the approximation accuracy of kinematic wave increases. According to Eq. (8) we observe that to obtain a numerical diffusion comparable with hydraulic one in diffusive wave, Dx ! ∞. However, from the theory of numerical methods it is known that the numerical grid allows to reproduce the waves longer than 2Dx only. Generally speaking, the different methods lead to a solution that represents the long wave properties only, whereas the short waves are considerably disturbed. Therefore if the flood wave has large gradients, the oscillations of the solution should be
R. Szymkiewicz / Journal of Hydrology 216 (1999) 248–251
expected even when Dx satisfies Eq. (8). For this reason it is rather impossible to obtain some general formulas defining the optimal values of distance step Dx and time step Dt for Muskingum model because they can depend on the character of flood wave also. These problems are connected with a general principle; the increasing order of approximation accuracy does not have to improve the accuracy of solution when the grid is sparse. Inversely, the solution accuracy suffers. The remarkable improvement of accuracy of the solution occurs when the increasing order of approximation appears together with decreasing mesh dimensions (Fletcher, 1991). Unfortunately it is not interesting for users of hydrological models because this process, reducing the numerical error
251
leads to the exact solution of kinematic wave which, as is known, for C const corresponds to pure translation of flood wave without any shape deformation.
References Abbott, M.B., Basco, D.R., 1989. Computational Fluid Dynamics. Longman Scientific and Technical. Cunge, J., 1969. On the subject of a flood propagation computation method (Muskingum method). J. of Hydraulic Research 7, 205– 230. Fletcher, C.A., 1991. Computational Techniques for Fluid Dynamics, 1, Springer-Verlag, Germany.