CIVE 530 - OPEN-CHANNEL HYDRAULICS LECTURE 12: FLOOD ROUTING

12.1  ROUTING OF KINEMATIC WAVES

The kinematic wave equation is obtained by combining the water continuity equation with a statement of uniform flow (the Manning or Chezy equations) (Ponce: Page 279) The water continuity equation is: ∂Q/∂x + ∂A/∂t = 0 The statement of uniform flow is: Sf - So = 0 Also: Q = α Aβ Taking the derivative: ∂Q/∂A = α β Aβ - 1 Therefore, the derivative is: ∂Q/∂A = β Q/A = β V = c in which c = kinematic wave celerity, or Seddon celerity (Seddon's Law), also expressed as follows: c = ∂Q/∂A = (1/T) ∂Q/∂y Multiplying this equation by the water continuity equation leads to the kinematic wave equation: ∂Q/∂t+ c ∂Q/∂x = 0 This equation can be solved numerically in a variety of ways, linear (constant c) or nonlinear (variable c). The linear centered-in-time/centered-in-space scheme (weighting factor 0.5 in space and time) is accurate, but unstable for any Courant number other than 1, where C is (Ponce: Page 281): C = c Δt / Δ x The linear backward-in-time/backward-in-space scheme is not accurate (it introduces too much numerical diffusion), but it is stable for any Courant number C. (Ponce: Page 283) The linear forward-in-time/backward-in-space scheme is not accurate (it solves the diffusion wave equation), and unstable for any Courant number C > 1. The linear backward-in-time/forward-in-space scheme is not accurate (it solves the diffusion wave equation instead), and unstable for any Courant number C < 1. These two schemes are used in the kinematic wave solution of HEC-1 and HEC-HMS. The SCS convex method of flood routing is a variant of the linear forward-in-time/backward-in-space scheme, where C is determined empirically. (Ponce: Page 283) There is no practical method to solve the kinematic wave equation accurately and without instability. A little off-centering of the derivatives (to weighting factor 0.51 to 0.55) would produce stability through a wide range of Courant numbers, but at a cost of the introduction of a small amount of numerical diffusion.

12.2  ROUTING OF DIFFUSION WAVES

The diffusion wave equation is obtained by combining the water continuity equation with a statement of conservation of momentum under gradually varied steady flow. The water continuity equation is: ∂Q/∂x + ∂A/∂t = 0 The statement of conservation of momentum for gradually varied steady flow : Sf = So - ∂d /∂ x Combining these two equations into one leads to the diffusion wave equation (Ponce: Page 290): ∂Q/∂t + c ∂Q/∂x = ν∂Q2/∂x2 in which ν = hydraulic diffusivity, or Hayami diffusivity, expressed as follows: ν = qo / (2 So) This equation is of second order; therefore, it can describe physical diffusion. Compare this equation with the kinematic wave equation, which is of first order and therefore, cannot describe physical diffusion. We conclude that the kinematic wave equation can only describe diffusion through the introduction of numerical diffusion. On the other hand, the diffusion wave equation can describe diffusion on its own, without the need to introduce numerical diffusion. The diffusion wave equation can be solved numerically using the implicit six-point Crank-Nicolson scheme, used for the numerical solution of second-order parabolic equations. The diffusion wave equation can also be solved numerically in an approximate form, by matching the physical diffusion of Hayami with the numerical diffusion of Cunge. This gives rise to the Muskingum-Cunge method.

12.3  MUSKINGUM-CUNGE METHOD

The Muskingum method of flood routing was developed in 1938 by McCarthy, who was working on the flood control schemes in the Muskingum river basin, in Ohio. The method is based on specifying a relationship in a reach between inflow, outflow and storage, such that (Ponce: Page 272): S = K [XI + (1-X)O] The parameter K is interpreted as the travel time through the reach, defined as follows: K = Δx/c The parameter X is interpreted as a weighting coefficient, limited in the range 0. ≤ X ≤ 0.5. In the Muskingum method, X is determined by estimation from similar data (hydrologically homogeneous watersheds), or empirically, by using gage data. Values of X = 0.2 and 0.3 have been used in practice. K controls the translation of the flood hydrograph; X controls the diffusion. In 1969, Cunge linked the value of X to the numerical diffusion of the problem at hand, defined in terms of Δx and Δt. The numerical diffusion coefficient derived by Cunge is (Ponce: Page 293): νn = c Δx [(1/2) - X] Matching physical (Hayami's) and numerical (Cunge's) diffusion coefficient leads to a predictive equation for X (Ponce: Page 294): X = (1/2) [1 - qo/(SocΔx)] Defining the cell Reynolds number D as follows (Ponce: Page 296): D = qo/(SocΔx) leads to the Muskingum-Cunge routing equation: Qj+1n+1 = C0Qjn+1 + C1Qjn + C2Qj+1n With the routing coefficients: C0 = (-1 + C + D) / (1 + C + D) C1 = (1 + C - D) / (1 + C + D) C2 = (1 - C + D) / (1 + C + D) Given c, Δx, Δt, qo, So, one can calculate C and D and proceed with the routing. c and qo have to be estimated either linearly (a value of qo and a corresponding value of c), or nonlinearly. Linear calculation: Ponce: Page 295. If the calculations are nonlinear, the routing parameters C and D are allowed to vary with the flow. Example: Route a flood wave with the following characteristics: peak flood Qp= 1000 m3/s; baseflow Qb= 0 m3/s; So = 0.000868; Ap = 400 m2; Tp = 100 m; β = 1.6; Δx = 14.4 km; Δt = 1 hour.   Solution: For lack of other data, base calculation on peak discharge. The mean velocity V = Qp/Ap = 1000/400 = 2.5 m/s. The wave celerity is: c = β Vp = 1.6 × 2.5 = 4 m/s. The flow per unit width qo (based on peak discharge) = Qp/Tp = 1000/100 = 10 m2/s. The Courant number C = c Δt/Δx = 4 m/s × 3600 s / 14400 m = 1. The cell Reynolds number D = qo/(SocΔx) = 10 / (0.000868 × 4. × 14400) = 0.2. The routing coefficients are: C0 = 0.091; C1 = 0.818; and C2 = 0.091. The routing calculations are shown below. Time Inflows C0I2 C1I1 C2O1 O 0 0 - - - 0 1 200 18.2 0.0 0.0 18.20 2 400 36.4 163.6 1.66 201.66 3 600 54.6 327.2 18.35 400.15 4 800 72.8 490.8 36.41 600.01 5 1000 91.0 654.4 54.60 800.00 6 800 72.8 818.0 72.80 963.60 7 600 54.6 654.4 87.69 796.69 8 400 36.4 490.8 72.50 599.70 9 200 18.2 327.2 54.57 399.97 10 0 0.0 163.6 36.4 200.00 11 0 0.0 0.0 18.2 18.2 12 0 0.0 0.0 1.66 1.66 13 0 0.0 0.0 0.16 0.16 Note that the peak outflow is 963.6 and it occurs at time 6 hr. The wave has diffused from a peak of 1000 to 963.6, and has translated from t = 5 hr to t = 6 hr.

12.4  DYNAMIC WAVE ROUTING

The Saint-Venant equations can be solved implicitly by using the Preissmann scheme. There is a need to specify the scheme's weighting factor, commonly set at 0.55 or 0.6. This weighting factor serves the purpose of introducing enough numerical diffusion so that the inherent nonlinearities, which could lead to instabilities, could be controlled (attenuated). The solution can be linearized, or solved in a nonlinear fashion using the Newton-Raphson method. The linear solution is used as first guess of the Newton-Raphson iteration. The Preissmann scheme results in a solution requiring the inversion of a diagonal matrix. The double-sweep algorithm is used to solve the matrix. The solution is not complete until the downstream boundary condition has been properly specified. Discharge is usually specified upstream (coming from the hydrology). Stage is usually specified downstream. A rating is used in lieu of downstream stage. The specification of a single-valued rating (kinematic), although practical and common, contradicts the solution at the boundary: the wave cannot be dynamic at the interior points and kinematic at the downstream boundary. A way out of this difficulty is to specify a looped rating, by taking into account the water surface slope at the downstream boundary. Another way is to artificially extend the downstream boundary a length sufficiently downstream (2 to 3 times the reach length) so that the error imposed by the kinematic downstream boundary condition does not affect the calculation at the actual (physical) reach's interior points.

050509