Question 3 of 7: Hydrologic modelling, reservoir and lake routing

Problem 3: Hydrologic modelling, reservoir and lake routing (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.

Part (i) — Two model equations, side by side (10 marks)

The question names two building blocks of every routing model. The first, the hydrologic transport model, is the lumped continuity or storage equation that transports a flood wave through a reach or a reservoir by mass conservation alone. The second, the non-linear outflow-storage equation, is the relation that must be supplied alongside continuity to close the problem, and which is non-linear whenever the outlet is a weir, an orifice or a natural channel control.

(1) Hydrologic transport model. The governing equation is the lumped continuity equation applied to the reach or reservoir taken as a single control volume,

which, integrated over a routing interval $\Delta t$ and written in the finite-difference form used for computation, becomes

where $I$ is inflow, $O$ outflow and $S$ storage, all in consistent units (m³/s and m³, with $\Delta t$ in seconds). This single equation contains two unknowns at the new time level, $O_2$ and $S_2$, so a second relation — the storage-outflow relation — is always required.

(2) Non-linear outflow-storage equation. The general form is

with $K$ a storage coefficient and $m$ an exponent. When $m = 1$ the element is a linear reservoir, $S = KO$, and $K$ has the dimension of time. When $m \neq 1$ the relation is non-linear, and this is the normal case for a real outlet, because the hydraulic control is itself non-linear: a sharp-crested or broad-crested weir gives

and a submerged orifice gives $O = C_d A_o\sqrt{2gH}$, that is, $O \propto H^{1/2}$, while the storage above the outlet in a pond of surface area $A_s$ is $S \approx A_s H$. Eliminating the stage $H$ between the two gives $S = K O^{2/3}$ for the weir and $S = K O^{2}$ for the orifice. For a natural lake or a river reach at a fixed control the same idea appears as $S = a Q^{b}$ fitted to measurements.

The two are used together. Rearranging continuity so that all unknowns stand on the left produces the storage-indication form,

in which the right-hand side is entirely known at each step. A curve of $\left(2S/\Delta t + O\right)$ against $O$ is prepared once from the stage-storage and stage-discharge relations, and the new outflow is read from it — which is exactly where the non-linear outflow-storage equation enters the calculation.

A short numerical illustration of the non-linearity: a detention pond of surface area 12 000 m² discharges over a 6.0 m broad-crested weir with $C = 1.7$. At a head of $H = 0.80$ m above the crest, the outflow is $O = 1.7 \times 6.0 \times 0.80^{3/2} = 7.30$ m³/s while the storage above the crest is $S = 12\,000 \times 0.80 = 9600$ m³. Doubling the storage to 19 200 m³ ($H = 1.60$ m) does not double the outflow: it raises it to $1.7 \times 6.0 \times 1.60^{3/2} = 20.6$ m³/s, a factor of 2.83, because $O \propto S^{3/2}$ here. That disproportion is precisely why a pond attenuates.

Part (ii) — Three key steps in the field verification of a rainfall-runoff model (5 marks)

Step 1 — Establish an independent field data set that was not used in calibration (split-sample design). Verification is meaningless if the data confirming the model are the data that shaped it. The record is therefore divided in advance: one period, containing a range of event magnitudes and antecedent conditions, is used to fit the parameters, and a separate period is reserved untouched for verification. In the field this means installing or confirming the instrumentation that will supply the verification data — a tipping-bucket rain gauge network dense enough to resolve the storm pattern, a stream gauge with a current rating (ideally a Water Survey of Canada station with published HYDAT data), and, where the model contains a soil-moisture or snow store, soil-moisture probes and snow-course measurements. Where the model is to be applied to a changed condition, such as post-development land use, verification should include at least one event observed under that condition or a hydrologically analogous gauged basin.

Step 2 — Run the model on the reserved events with the calibrated parameters frozen, and compare simulated with observed hydrographs on multiple criteria. No parameter may be touched at this stage; that is the whole point. The comparison must go beyond a single summary statistic, because a model can match volume while getting timing badly wrong. In practice the criteria are: total runoff volume (percent error); peak discharge (percent error); time to peak (hours of lag); the shape of the recession; and an integral goodness-of-fit measure such as the Nash–Sutcliffe efficiency or the root-mean-square error. Plotting simulated against observed hydrographs event by event, and simulated against observed peaks across all events, exposes systematic bias — a model that consistently under-predicts large events and over-predicts small ones has a loss-model problem rather than a routing problem.

Step 3 — Diagnose residual error against field evidence, and confirm the internal states, not just the outlet discharge. A model can reproduce the outlet hydrograph for the wrong reasons — too much infiltration compensated by too little evapotranspiration, for example. The third step is therefore to check the model’s internal variables against whatever field observations exist: measured soil-moisture profiles against the simulated soil store, observed water levels in the stormwater pond against the routed stage, high-water marks and flooded extents from a post-event field walk against the simulated peak stage, observed snow-water equivalent against the simulated snowpack, and baseflow separation against measured groundwater levels. Any systematic discrepancy is traced to a specific process, and if the model must be revised, the revision requires a fresh verification on data not used in the revision. The outcome of this step is a stated range of applicability — the event magnitudes, seasons and land-use conditions for which the model has been shown to work — which is what a professional engineer signs.

Part (iii) — Fundamentals of the Muskingum routing method (5 marks)

The Muskingum method, developed for the Muskingum River flood-control studies, is a hydrologic (storage) routing method for a river reach, and it is the standard answer to the closely related reservoir problem because it reduces to level-pool routing when the wedge term vanishes. Its fundamental idea is that the storage in a river reach cannot be described by outflow alone, as it can in a level-pool reservoir, because during a flood the water surface in the reach is not horizontal.

Storage is therefore split into two parts. Prism storage is the volume that would exist under a water surface parallel to the bed at the steady discharge corresponding to the outflow, and is taken proportional to outflow, $K O$. Wedge storage is the additional volume held between that parallel surface and the actual sloping surface; on the rising limb, inflow exceeds outflow, the surface tilts up toward the upstream end and the wedge is positive, while on the falling limb the wedge is negative. It is taken proportional to the difference between inflow and outflow, $K X (I - O)$. The total is

Here $K$ is the storage-time constant, approximately the travel time of the flood wave through the reach (dimension T), and $X$ is a dimensionless weighting factor between 0 and 0.5 that measures how much the inflow contributes to storage: $X = 0$ makes storage a function of outflow alone, which is the level-pool reservoir case, while $X = 0.5$ gives equal weight to inflow and outflow and yields pure translation with no attenuation. For natural channels $X$ typically lies between 0.1 and 0.3, with 0.2 a common first estimate.

Substituting this storage relation into the finite-difference continuity equation and solving for the unknown outflow gives the working equation,

with the routing coefficients

and the identity $C_0 + C_1 + C_2 = 1$, which must be checked at every application because it guarantees that mass is conserved. The parameters are determined either from observed inflow and outflow hydrographs — plotting accumulated storage against the weighted discharge $XI + (1-X)O$ for trial values of $X$ and selecting the $X$ that collapses the loop into the narrowest single line, whose slope is then $K$ — or, for an ungauged reach, from the Muskingum–Cunge variant in which $K = L/c$ with $c$ the kinematic wave celerity and $X = \tfrac{1}{2}\left(1 - q/(S_0 c L)\right)$ computed from channel geometry and Manning roughness.

Two stability conditions must be respected: the routing interval must satisfy $2KX \le \Delta t \le 2K(1-X)$, and $\Delta t$ should be no greater than about one-fifth of the time to peak of the inflow hydrograph. If $\Delta t < 2KX$ then $C_0$ becomes negative, and the computed outflow can dip below the initial flow at the start of the event — the well-known negative initial outflow that signals a badly chosen interval.

Given. Reach parameters $K = 2.3$ h and $X = 0.20$ (assumed, typical of a natural channel), routing interval $\Delta t = 1.0$ h, inflows $I_1 = 50$ and $I_2 = 80$ m³/s, initial outflow $O_1 = 50$ m³/s.

Find. The routing coefficients, the stability check and the outflow at the end of the first interval.

1. Check the stability limits on the routing interval. With $2KX = 2(2.3)(0.20) = 0.92$ h and $2K(1-X) = 2(2.3)(0.80) = 3.68$ h, the requirement $0.92 \le \Delta t \le 3.68$ is satisfied by $\Delta t = 1.0$ h, so $C_0$ will be positive.
2. Evaluate the routing coefficients. The common denominator is $2K(1-X)+\Delta t = 3.68 + 1.0 = 4.68$ h, so $$\begin{aligned} C_0 &= \frac{1.0 - 0.92}{4.68} = 0.0171 \\ C_1 &= \frac{1.0 + 0.92}{4.68} = 0.4103 \\ C_2 &= \frac{3.68 - 1.0}{4.68} = 0.5726 \end{aligned}$$ Their sum is 0.0171 + 0.4103 + 0.5726 = 1.0000, so mass will be conserved.
3. Route one interval. Substituting into the working equation, $$O_2 = 0.0171(80) + 0.4103(50) + 0.5726(50) = \boxed{50.5\ \text{m}^3/\text{s}}$$ The outflow has barely begun to rise while the inflow has risen by 30 m³/s, which is the attenuation and lag the method exists to reproduce; repeating the step with each successive pair of inflows generates the whole outflow hydrograph.