War. Res. Vol. 26, No. 12, pp. 1563-1570. 1992 0043-1354/92 $5.00+0.00 Printed in Great Britain. All rights reserved Copyright .~ 1992 Pergamon Press Ltd IMPROVING THE NUMERICAL MODELING OF RIVER WATER QUALITY BY USING HIGH ORDER DIFFERENCE SCHEMES A. I. STAMOU @ Department of Civil Engineering. National Technical University of Athens, lroon Polytechniou 5, 15773 Zografou, Athens, Greece (First receired Not'ember 1991; accepted in revised form May 1992) Abstract--A one-dimensional numerical model is presented for the description of water quality in rivers. The model is an improvement over existing models, since it involves two differential schemes of third-order truncation error (QUICK and QUICKEST). which avoid the stability problems of central difference, while remaining relatively free of the inaccuracies of numerical dispersion associated with the upwind difference scheme. The model is successfully applied to four simple, steady-state and transient water quality problems. Results show that the commonly used upwind scheme performs very poorly. The central difference scheme performs satisfactorily in steady-state implicit calculations, but not in transient problems, where it shows numerical dispersion and significant oscillations. Calculations with QUICK and QUICKEST are stable and accurate involving the minimum error and minimum level of numerical dispersion in both steady-state and transient problems. Key words--water quality, river modeling, difference schemes. QUICK. QUICKEST, numerical dispersion INTRODUCTION The majority of existing river water quality models are mechanistic and are based on differential equations derived by applying mass balances (see O'Connor et al., 1973) for the water quality variables of interest (e.g. BOD, DO, N, P). These equations have analytical solutions only for specific boundary conditions, as in the cases of the deterministic model of Miller and Jennings (1979) and the stochastic model of Thayer and Krutschkoff (1967). In most of the practical cases, however, an analytical solution does not exist and it is thus necessary to solve the differential equations numerically using the finite difference technique (Elliot and James, 1984a, b). Essentially, the finite difference technique attempts to approximate the continuous solution derived by the exact equations at a discrete number of points in space (grid points), which overlay the physical do- main considered (Fig. I) and in time (time levels), by the so-called difference equations. Difference equations are derived from the differential equations, when the latter are integrated over finite difference control volumes, defined by the grid points (see Fig. 1). In the integration procedure, certain difference schemes are applied for the interpolation of variables between grid points. In the majority of existing water quality models encountered in the literature ('l'homann and Mueller, 1987; Elliot and James, 1984a, b) the simple low-order upwind (UDS) and central difference (CDS) schemes are applied, which assume a zeroth and first-order (linear) variation, respectively. Despite their simplicity, these schemes have signifi- cant disadvantages. The UDS, although numerically stable, involves a first-order truncation error, which introduces a significant amount of numerical dis- persion and results in severe inaccuracies. The CDS is more accurate than the upwind procedure; since it involves a second-order truncation error, but it may become unstable in steady-state calculations and show numerical dispersion in transient calculations. To overcome these problems several investigators have proposed various numerical schemes, which are mainly of first or second order, in the sense that they involve one or two grid points in the interpolation procedure. For example, Bella and Grenney (1970) used the UDS, in which they introduced a pseudo- dispersion to reduce the numerical dispersion.Holly and Preissmann (1977) presented their "two-point higher-order method" and concluded that this scheme is superiorto the multiple-pointhigher-order schemes. In the present work two quadratic interp- olation numerical schemes are proposed, which are based on the QUICK (quadratic upwind interp- olation for convection kinematics, QDS) and QUICKEST (quadratic upwind interpolation for convection kinematics with estimated streaming terms, QEDS) schemes introduced by Leonard (1979). These schemes, which involve a third-order truncation error, are more accurate than CDS and they also avoid its stability problem. Furthermore, they remain relatively free of the inaccuracies of numerical dispersion which are associated mainly with the UDS. 1563