USGS -- SMIG --
Surface-water quality and flow Modeling Interest Group

Simulation Model for Open-Channel Flow and Transport

by Raymond W. Schaffranek

USGS, Water Resources Division
12201 Sunrise Valley Drive, MS 430
Reston, VA 20192
Internet: rws@usgs.gov
Phone: (703) 648-5891
FAX: (703) 648-5484

Editor's note:
This paper was written for the First Federal Interagency Hydrologic Modeling Conference, held in Las Vegas on April 19-23, 1998.

Citation:
Schaffranek, R.W., 1998, Simulation Model for Open-Channel Flow and Transport, in Proceedings of the First Federal Interagency Hydrologic Modeling Conference, April 19-23, 1998, Las Vegas, NV: Subcommittee on Hydrology of the Interagency Advisory Committee on Water Data, p. 1-103 to 1-110.

Contents

Abstract

A flow-simulation model, formulated using an extended form of the one-dimensional de Saint Venant equations of unsteady open-channel flow, has been augmented to include solution of the advection-dispersion transport equation. The weighted, four-point implicit, finite-difference approximation of the unsteady-flow equations permits solutions at large time steps. The model is fully capable of simulating unsteady flow throughout a network of open channels connected in a dendritic or looped pattern. The flow model accommodates dynamic tributary flows and controlled diversions as well as lateral inflows and overbank storage. A mixed Eulerian/Lagrangian approach is used to solve the one-dimensional advection-dispersion equation for coupled simulation of solute transport. The transport solution technique avoids numerical dispersion associated with conventional finite-difference methods for advection-dominated problems. The technique is free of Courant number restrictions and can be applied using large time steps consistent with the flow-model discretization. Coupled solution of the unsteady flow and transport equations circumvents interpolation errors inherent in de-coupled simulation approaches and enables concurrent treatment of longitudinal pressure differentials due to variable density gradients.

Introduction

A Lagrangian particle-tracking scheme, previously developed and incorporated in the branch-network dynamic flow model, BRANCH, (Schaffranek et al. 1981, Schaffranek 1987a) permits examination of the retention times of parcels of water in open-channel networks. This particle-tracking scheme was first demonstrated in application to the tidal Potomac River (Schaffranek 1987b). The model was used to investigate the flushing capacity and retention properties of the tidal-river system for purposes of analyzing factors contributing to the development of algal blooms and the fate of phytoplankton. Retention times were found to vary considerably in response to changes in freshwater inflows, tidal dynamics, and meteorological conditions. Use of the model made it possible to identify local flow patterns under various combinations of boundary conditions and to investigate the role of tidal trapping in the fate of phytoplankton. Retention times in the main Potomac River channel were found to vary from weeks to months for moderate to low inflows, but to be of only several days duration for high inflows produced by upstream storm events.

In the effort described herein, the Lagrangian particle-tracking scheme of the dynamic flow model has been extended to include full solution of the transport equation by incorporating treatment of the physical process of dispersion. This direct coupling of the unsteady-flow and transport solutions eliminates the errors of numerical dispersion inherently introduced in de-coupled approaches wherein flow velocities must be interpolated to comply with the grid-point requirements of the transport scheme. Moreover, a coupled flow/transport simulation approach is required in situations where the fluid density varies longitudinally---as typically encountered in tidal-influenced, sediment-laden, or thermally-varying riverine flows---necessitating treatment of pressure differentials in the unsteady-flow equations.

In this paper, the BRANCH model of the tidal Potomac River network (Schaffranek 1987b) is used to demonstrate the coupled flow/transport simulation capability of the generic model. A previous numerical simulation, based on tracking parcels of water representing a conservative, non-dispersive constituent (Schaffranek 1987b), is repeated and revised to illustrate the capabilities of the model in simulating the combined effects of advection and dispersion on solute transport. Numerical simulations, conducted at varied Peclet numbers, illustrate the stability of the method when advection is the dominant process and the controlled response of the method for representation of physical dispersion without the introduction of numerical artifacts.

Unsteady Flow Equations

Unsteady flow in open channels is governed principally by dynamic equilibrium of momentum changes due to inertial and convective accelerations, differential pressure forces, and gravitational and shear-stress effects of the bed and friction slopes. A variety of differential-equation sets can be derived for depicting unsteady open-channel flow (Lai 1986). The equation set employed in this model uses flow discharge, Q, and water-surface elevation above a horizontal datum (stage), Z, as the pair of dependent variables and longitudinal distance along the channel, x, and time, t, as the two independent variables:

equation      (1)

equation      (2)

In derivation of these equations, a local Cartesian coordinate system was used in which x, y, and z are the streamwise, transverse, and vertical axes, respectively, and the flow discharge is normal to the y-z plane in the positive x direction. Assumptions used in their formulation are that the water is of homogeneous density, the vertical pressure gradient is hydrostatic, and the channel bottom is rigid---or relatively stable and fixed with respect to time---with a mild uniform slope. (The longitudinal pressure gradient term is not included in the equation set for this analysis.) Other symbols in the above equations, B, q, beta, A, g, Sf, u prime, xi, and Vw, represent the channel top width, lateral inflow per unit length of channel, the Boussinesq momentum-correction coefficient, the cross-sectional area, gravitational acceleration, the friction slope, the x-component of lateral-inflow velocity, the wind-stress coefficient, and the wind velocity (occurring at an angle alpha to the longitudinal axis of the channel), respectively. The momentum coefficient, beta, defined as glyph u2dA/U2A, in which u is the velocity of water flowing through some finite elemental area dA and U is the mean flow velocity in the cross-sectional area A, is used to adjust for any nonuniform velocity distribution in the channel cross section. A Manning formulation of the friction-slope term is used in which Sf takes the form (eta/lambda)2Q |Q| /A2R4/3, wherein lambda = 1.0 for metric or 1.486 for foot-pound units, R is hydraulic radius, and the symbol eta is used in place of the Manning n to indicate that the coefficient is being used to represent frictional resistance under unsteady-flow conditions. The dimensionless wind-stress coefficient, eta, is a function of the water-surface drag coefficient, Cd, the water density, rho, and atmospheric density, rho-a, expressed as Cd(rho-a/rho).

Four-Point Implicit Solution

The BRANCH model (Schaffranek et al. 1981, Schaffranek 1987a) uses a four-point, implicit, finite-difference approximation of the unsteady-flow equations (1) and (2). In the four-point, Preissmann (1961) or box scheme, spatial derivatives are centered in space and weighted in time according to a weighting factor theta. The finite-difference approximation of the spatial derivatives can vary from box-centered (theta = 0.5) to fully forward (theta = 1), which is the range for which its stability has been proven. See, for example, Fread (1974) for stability analysis using a linear wave approximation with a linearized form of the friction-slope term, and Samuels and Skeels (1990) for Fourier stability analysis of the linearized numerical equations including the general form of the friction term. To avoid computational-mode oscillations, a practical lower limit for theta appears to be 0.6 (Schaffranek et al. 1981). Time derivatives are centered in space; geometric properties and functional quantities are centered in space and time weighted according to theta. After the flow equations are coupled recursively to eliminate internal nodes, the matrix of flow and boundary-condition equations is solved in iterative fashion by Gaussian elimination using maximum pivot strategy.

Transport Simulation Component

The model solves the one-dimensional, advection-dispersion, transport equation describing the concentration of a solute as a function of time and distance, i.e.,

equation      (3)

in which C(x,t) is the cross-sectional average concentration of a solute at location x along the longitudinal axis of channel at time t, U is the mean cross-sectional flow velocity (Q/A), and Dx is the longitudinal dispersion coefficient. Solution of the transport equation is directly coupled to the unsteady-flow model.

Transport Solution

The unknown flow conditions Q(x,t) and A(x,t) are provided directly by solution of the unsteady-flow equations (1) and (2); however, the transport equation (3) is not solved entirely by a Eulerian finite-difference technique. Instead, a mixed Eulerian/Lagrangian approach is used in which the advection and dispersion processes are treated separately and distinctly. The more difficult advection process (u dc/dx) is resolved by tracking particles representing parcels of water in a Lagrangian reference frame, whereas the dispersion process (d dc/dx) is solved in a Eulerian framework by a finite-difference technique. A cubic spline is employed to interpolate the concentration profile of advected particles. An explicit, space-centered finite-difference scheme is used to represent the dispersion process. This solution approach avoids the well-known difficulties associated with numerical dispersion in advection-dominated systems introduced by pure Eulerian finite-difference solution methods. (See, for example, Gray and Pinder (1976).)

The cubic spline interpolant for the concentration profile is based on the method developed by Akima (1978). This spline interpolation technique has the desirable feature of minimizing the development of oscillations manifested as over- and under-shoot in the interpolated function. The Akima spline attempts to preserve the shape of the concentration profile as reflected by the data. The behavior of the Akima spline method is compared with a natural cubic spline (de Boor 1978) interpolation of the same function, potentially representative of a temporal or spatial constituent concentration profile, in Fig. 1. Endpoint conditions are automatically determined by the Akima method, whereas values of zero are specified for second derivatives at the endpoints for the natural spline. The extraneous inflection points of the natural cubic spline (Fig. 1a), yielding over- and under-shoot in the interpolated function, are clearly illustrated in contrast to the smoother interpolant produced by the Akima spline method (Fig. 1b). Moreover, the Akima spline reduces the potential for the generation of negative concentrations as is evidenced in the interpolation produced by the natural spline (Fig. 1a). These attributes of the Akima spline interpolation method contribute to improved conservation properties for the transport-solution technique.

figure 1a (fig. 1a)

figure 1b (fig. 1b)

Figure 1. Natural (a) and Akima (b) spline interpolants of hypothetical concentration profile.

The Tidal Potomac River

The tidal-river segment of the Potomac River extends downstream from the head-of-tide, near Chain Bridge in the northwest quadrant of the District of Columbia, to Indian Head, Maryland, a distance of nearly 50 km. The cross-sectional area of the tidal river expands more than fortyfold between Chain Bridge and Indian Head, increasing from 232 to 9,960 m2. The corresponding width increases from 44 to 1,950 m. Although the channel bottom is somewhat irregular depths, in general, range from about 9 m at Chain Bridge to 12 m at Indian Head. Several tributaries and tidal inlets, in which depths are typically 3 m or less, adjoin the main channel to form the tidal-river network.

Model Implementation

The model of the tidal Potomac River system was originally developed (Schaffranek 1987b) in order to study the flow dynamics of the network of channels. The model implementation consists of 25 branches joining or terminating at 25 nodes. A total of 66 cross sections, at intervals ranging from 0.7 to 4.8 km, are used to depict the irregular geometry of the tidal river and its tributaries and tidal inlets.

Boundary conditions used to conduct numerical simulations include zero discharges at the ends of tributaries and inlets, tide elevations recorded at Indian Head, and freshwater inflows derived from a rated gaging station near the head-of-tide. Freshwater inflows at the head-of-tide average 323 m3/s, but have ranged from 3 to 13,700 m3/s. Tidal amplitude at Indian Head is about 55 cm, but wind effects can significantly dampen tidal propagation. Wind speed and direction data collected at Indian Head are used to evaluate wind-forcing effects. Flow simulations are conducted using a 15-minute time step, which yields a Courant number Cr ~ 6. The model was calibrated and verified using five sets of tidal-cycle flow discharges measured at three different locations within the tidal-river system. Measured and model-computed flood and ebb flow volumes agree within 10 percent (Schaffranek 1987b).

Numerical Simulation Results

Flow data from the Potomac River, for the 30-day period beginning August 15 and ending September 13, 1981, are used to demonstrate the coupled flow and transport simulation capabilities of the model. These data are of particular interest due to the extremely low freshwater inflows that prevailed during the period. Inflows were less than the 5-year (1979-1983) average and less than 100 m3/s, except for the last six days of the period. Ebb and flood tidal cycle discharges at Indian Head, as computed by the flow model, ranged from approximately 4,400 to 4,700 m3/s.

Particle Tracking: The transport of a conservative-type substance by advective processes alone is illustrated by the simulation results presented in Fig. 2. These results were obtained by tracking index particles representing parcels of water. In the simulation, nine index particles---potentially representing the location of a conservative, neutrally buoyant, passive solute or suspended particulate matter---were tracked through the main channel of the tidal river. The paths of travel of these particles are shown in the central graph of the figure. The vertical axis of the time-of-travel graph represents the main tidal Potomac River channel between Chain Bridge and Indian Head. The upper hydrograph in the figure shows the freshwater inflow recorded near Chain Bridge and lower hydrograph presents the tidal water-surface elevations simultaneously recorded at Indian Head. Together these hydrographs constitute the instantaneous boundary conditions used to solve the unsteady-flow equations.

(fig. 2)

Figure 2. Track of injected particles due to pure advection in the Potomac River for August 15 through September 13, 1981. [A larger figure is available (64k GIF).]

The simulation results plotted in Fig. 2 demonstrate the combined effects of freshwater inflows and tidal forcing on the displacement of index particles and give indications of the retention properties and flushing characteristics of the tidal-river system. Throughout the period, the movement of index particles is variously influenced by the predominately semi-diurnal tides at Indian Head as they cycle through spring-neap variations and changing freshwater inflows at the head-of-tide. Near the end of the period, increased freshwater inflows, declining tidal elevations, and a significant wind event act to produce a pronounced downstream particle displacement. As is evident from the particle tracks, there is a tendency throughout the period for all particles to move closer to one another as they are advected along the channel. This variability in particle displacement can be attributed to the dynamic flow conditions and, in particular, to the magnitude and resultant dominance of the net tidal influx into the system as compared with the magnitude of freshwater inflows. This tidal dominance yields extended periods of little or no net downstream, and occasionally upstream, displacement of particles, which contributes to a diminished flushing capacity.

Constituent Transport: The transport of a conservative, non-reactive constituent along the tidal Potomac River, under the influences of both advection and longitudinal dispersion, is illustrated in Fig. 3. The results of four transport simulations are presented in Fig. 3b for the initial, hypothetical, constituent concentration profile presented in Fig. 3a. The transport simulation is for the same conditions as used to generate the particle-tracking results presented in Fig. 2. Distances on the horizontal axis of the figure are referenced to Chain Bridge. The four concentration profiles of a transported constituent, illustrated in Fig. 3b, present simulation results midway through the 30-day period at midnight on August 29. Values of 0, 1, 5, and 10 m2/s were assigned to the dispersion coefficient, Dx. These coefficient values represent pure advection, Dx = 0, and mixed advection-dispersion conditions, Dx = 1, 5, 10, identified by Peclet numbers, Pe = (peclet number), of 20, 4, and 2, respectively, for the mean flow velocity and 800 m grid-point spacing, deltax, of particles.

figure 3a (fig. 3a)

figure 3b (fig. 3b)

Figure 3. Initial constituent concentration (a) and transported concentrations after 15 days (b) for varied longitudinal dispersion coefficients, Dx = 0, 1, 5, 10.

As Fig. 3 illustrates, the transport solution of the model yields stability in all four numerical simulations. In particular, no instabilities are visible and no numerical dispersion is evident in the simulation results in which pure advection is considered. The simulations, in which dispersive transport is also evaluated, i.e., Dx = 1, 5, 10 m2/s, illustrate that the attendant spreading of the concentration profiles can contribute to a lengthening of the retention time for a constituent within the system. Thus, the compression of the concentration profile as demonstrated in the pure advection results of the particle-tracking simulation will be diminished, and potentially overwhelmed, by the process of physical dispersion with resultant consequences on the evaluation of transport of a constituent through the system. By full consideration of advection and dispersion transport processes such as this, a more comprehensive analysis of the flushing capacity and retention properties of the tidal-river network can be conducted. This will lead to more definitive evaluations of the fate and resultant effects of constituents, potentially introduced or naturally occurring, within the system.

Summary and Conclusions

The branch-network dynamic flow model, BRANCH, has been extended to include coupled solution of the advection-dispersion equation for simulating solute transport within an open-channel system. The mixed Eulerian/Lagrangian, transport simulation approach yields stable solutions that are free of numerical dispersion introduced in pure Eulerian methods. Simulations, using the previously calibrated model as implemented for the tidal Potomac River, demonstrate the capabilities of the flow/transport model in evaluating constituent transport in a network of open channels. This transport simulation capability can be used to investigate the fate of constituents, such as pollutants or contaminants, potentially introduced into an open-channel network. Coupled solution of the unsteady-flow and solute-transport equations permits evaluation and consideration of density effects on the flow and eliminates interpolation errors inherent in de-coupled flow/transport simulation approaches. The simulation model also can be used to investigate the retention time and flushing properties of coastal network systems composed of interconnected channels that are variously influenced by freshwater inflows, tidal forces, and the stochastic effects of weather fronts.

References

Akima, H. (1978). A method of bivariate interpolation and smooth surface fitting for irregularly distributed data points: ACM Transactions on Mathematical Software, 4, 148-159.

de Boor, C. (1978). A practical Guide to Splines: Springer-Verlag, New York, N.Y., 392 p.

Fread, D.L. (1974). Numerical properties of implicit four-point finite difference equations of unsteady flow: Technical Memorandum NWS HYDRO-18, U. S. National Oceanic and Atmospheric Admin., Silver Spring, Md., 38 p.

Gray, W.G., and Pinder, G.F. (1976). An analysis of the numerical solution of the transport equation: Water Resources Research, 12(3), 547-555.

Lai, C. (1986). Numerical modeling of unsteady open-channel flow: Advances in Hydroscience, V.T. Chow and B.C. Yen, eds., Vol. 14, Academic Press, Orlando, Fla., 163-333.

Preissmann, A. (1961). Propagation of translatory waves in channels and rivers; (in French), Proc., First Congress of French Association for Computation (AFCAL), Grenoble, France, 433-442.

Samuels, P.G., and Skeels, C.P. (1990). Stability limits for Preissmann's scheme: ASCE Journal of Hydraulic Engineering, 116(8), 997-1012.

Schaffranek, R.W. (1987a). Flow model for open-channel reach or network: Professional Paper 1384, U.S. Geological Survey, Denver, Colo., 11 p.

Schaffranek, R.W. (1987b). A flow-simulation model of the tidal Potomac River. Water-Supply Paper 2234-D, U.S. Geological Survey, Denver, Colo., 41 p.

Schaffranek, R.W., Baltzer, R.A., and Goldberg, D.E. (1981). A model for simulation of flow in singular and interconnected channels: Book 7, Chap. C3, Techniques of Water Resources Investigations, U.S. Geological Survey, Denver, Colo., 110 p.

Back to the SMIG Features Page

button bar

SMIG Home Mailing List Features Conferences Classes Reading Room Model Archives Feedback Home | Mailing List | Features | Conferences | Classes | Reading | Model Archives | Feedback

Stewart Rounds, SMIG coordinator <sarounds@usgs.gov>
U.S. Geological Survey
http://smig.usgs.gov/SMIG/features_0998/branch.html
Last modified Thursday, 10-Mar-2005 18:30:50 EST
Privacy Statement · Disclaimer · FOIA · Accessibility