Question 4 of 7: Channel and River Routing, and Flood-Wave Behaviour

Question 4: Channel and River Routing, and Flood-Wave Behaviour (20 marks)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

(i) Four concepts in applying the storage–budget concept to river routing (10 marks)

The question asks for one concept applied to one routing type, so the answer below develops the Storage-Budget (hydrologic) concept applied to river channel routing — that is, the Muskingum family of methods, in which the reach is treated as a storage element and the routing is done by accounting rather than by solving the equations of motion.

Concept 1 — the reach is a storage element, and continuity is the governing statement. The starting point is conservation of mass written for the volume of water held between an upstream and a downstream section. Over a routing interval Δt,

where I and O are the inflow and outflow discharges and S is the storage in the reach. Nothing here requires momentum: the method is a budget. The immediate physical consequences are the two features visible on Figure 4.1 — the outflow peak lags the inflow peak, and it is attenuated below it, because the reach fills while I > O and drains while O > I. The areas under the two hydrographs must be equal if there is no lateral inflow or loss, and checking that equality is the standard verification of a routing computation.

Concept 2 — a storage relation must be supplied to close the system. The continuity equation contains three unknowns (S₂, O₂ and, through them, the stage) but only one equation, so a second, independent relation between storage and discharge is required. For a reservoir this is easy: storage depends on stage alone, and stage determines outflow through the spillway rating, so S = f(O) is single-valued and level-pool routing applies. In a river reach it is not, because the water surface is sloped rather than horizontal: at a given outflow the reach holds more water on the rising limb than on the falling limb. The Muskingum relation captures this by making storage a weighted function of both ends of the reach,

with K the storage-time constant (approximately the travel time of the flood wave through the reach) and X a dimensionless weighting factor between 0 and 0.5. The term KO is the prism storage held beneath a surface parallel to the bed; the term KX(I − O) is the wedge storage created by the sloping surface, positive while the wave is rising and negative while it is falling. Setting X = 0 recovers the level-pool (reservoir) case.

Concept 3 — the parameters must be estimated from data, not assumed. Substituting the storage relation into continuity and solving for the unknown outflow gives the working equation

where the three coefficients must sum to unity — the arithmetic check that the routing conserves volume. When measured inflow and outflow hydrographs exist for a past flood, K and X are found graphically: accumulated storage is computed from continuity and plotted against the weighted discharge [XI + (1−X)O] for trial values of X; the value that collapses the loop into the narrowest single-valued line is adopted, and K is the slope of that line. Where no gauged pair exists, K is estimated as the reach length divided by the flood-wave celerity and X from reach hydraulics (typically 0.2–0.3 for natural channels, near 0 for a reach that behaves like a pond).

Concept 4 — discretise, march forward in time, and respect the stability limits. With the coefficients fixed, the routing proceeds step by step from a known initial outflow, each interval using the current and previous inflows and the previous outflow. The time step is not free: it must be short enough to define the inflow hydrograph — a common rule is Δt no greater than T_p/5 — and it must satisfy 2KX ≤ Δt ≤ 2K(1 − X) if negative outflows and oscillation are to be avoided. A long river is subdivided into sub-reaches so that the travel time of each is comparable to Δt, with the outflow of one becoming the inflow of the next and any tributary added at the junction. The final step is verification: confirm ∑C = 1, confirm equality of the inflow and outflow volumes, and confirm that the computed peak falls where the physics says it should — later and lower than the inflow peak.

(ii) Defining and solving the spring-thaw flood-wave problem (10 marks)

The problem is to predict, at every section downstream of a snowmelt-generated flood wave, the discharge Q(x,t), the mean velocity V(x,t) and the depth or wave height y(x,t). Because three unknowns are wanted and they vary in both space and time, the storage-budget approach of part (i) is no longer sufficient — it returns discharge only, at the reach outlet only. The problem is set up as a one-dimensional unsteady open-channel flow problem and solved as follows.

Define the system and generate the upstream boundary condition. Delineate the contributing basin above the reach of interest and establish the snowmelt input. The melt rate is computed either from the energy budget of the snowpack or, for operational work, from a degree-day relation M = C_m(T_a − T_b) calibrated for the basin, applied to the snow water equivalent measured on snow courses and by snow pillows and, in Canada, supplemented by provincial snow-survey bulletins. The melt depth is combined with any concurrent rainfall and reduced by the infiltration the frozen or saturated soil can accept — typically very little in a spring thaw, which is precisely why these events are severe — and the resulting excess is transformed into an inflow hydrograph at the head of the reach by a unit hydrograph or a runoff model. That hydrograph I(t) is the upstream boundary condition.

Assemble the channel data and the remaining boundary and initial conditions. Survey or extract from LiDAR a series of cross-sections along the reach, together with bed elevations, so that the geometry functions A(y), B(y) and R(y) and the bed slope S₀ are known at each section; assign Manning roughness values to the main channel and the overbank, calibrated against a past gauged flood if one is available. Set the downstream boundary condition — a known rating curve, a critical-depth control at a weir or a step in the bed, or a stage hydrograph if the reach discharges into a lake, reservoir or estuary. Set the initial condition as the steady pre-melt flow profile computed by a backwater (gradually varied flow) calculation. Add lateral inflow q_l for tributaries and for melt entering along the reach.

Write the governing equations. The dynamic-wave (Saint-Venant) system expresses conservation of mass and of momentum for gradually varied unsteady flow in one dimension:

with the friction slope supplied by Manning’s equation, S_f = n²Q|Q|/(A²R^4/3). These two equations in the two dependent variables Q and y close the problem; the velocity follows as V = Q/A and the wave height as the difference between the routed depth and the initial depth.

Choose the level of approximation deliberately. The momentum equation contains, from left to right, the local acceleration, the convective acceleration, the pressure (depth-gradient) term and the gravity–friction balance. Dropping the two acceleration terms gives the diffusion wave, which retains attenuation and is adequate for most river flood forecasting; dropping the depth gradient as well gives the kinematic wave, in which the wave translates without attenuation at celerity c_k = dQ/dA ≈ (5/3)V for a wide channel with Manning friction. The full dynamic wave must be retained where the reach is flat, where downstream backwater controls the flow — an ice jam, a reservoir, a tidal estuary — or where the wave is abrupt. A spring thaw on a Canadian river usually demands the full form, because ice-jam formation and release is exactly the backwater-dominated case.

Solve numerically and verify. Discretise the reach into computational sections and solve the coupled equations by an implicit finite-difference scheme — the four-point Preissmann box scheme is the standard, and is what HEC-RAS and MIKE 11 implement — marching forward in time from the initial profile. Explicit schemes require the Courant condition Δt ≤ Δx/(V + √gy); implicit schemes are unconditionally stable but still need a step small enough for accuracy. The output is Q, y and V at every section and every time step, from which the peak stage, the arrival time and the flood-plain extent follow directly.

Verification and the physical checks matter as much as the solution. Confirm that volume is conserved between the upstream and downstream boundaries, that the outflow peak is later and lower than the inflow peak, and that the computed stage–discharge relation at a section forms the expected counter-clockwise loop rather than a single-valued curve — the loop is the direct signature of the unsteady terms and shows that the wave is being represented rather than a sequence of steady states. Finally, calibrate against a past thaw event if a gauged record exists, and carry a sensitivity analysis on Manning n and on the melt rate, since those two are always the least certain inputs.