Question 4 of 7: Hydrologic Modelling, Model Validation and Reservoir Routing

Question 4: Hydrologic Modelling, Model Validation and Reservoir 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.

(i) The water balance equation as a basic hydrologic model (8 marks)

The water balance is the statement of conservation of mass for a defined control volume over a defined period: inflow minus outflow equals the change in storage. Applied to a watershed bounded by its drainage divide, with the soil column and shallow aquifer included, it reads

$$P - R - G - E - T = \Delta S$$

with all terms expressed as depths over the basin area (mm) for the same interval: $P$ precipitation, $R$ surface runoff leaving at the outlet, $G$ net subsurface outflow across the boundary, $E$ evaporation, $T$ transpiration and $\Delta S$ the change in storage (soil moisture, snowpack, groundwater, surface water). Written for a general control volume it is simply

$$\frac{\mathrm{d}S}{\mathrm{d}t} = I(t) - O(t)$$

and every hydrologic model — from a one-line annual balance to a distributed continuous-simulation code — is this equation applied to a chosen set of stores at a chosen time step. That is what makes it the basic model: choosing the control volume, the time step and the number of stores is the modelling decision.

How it is used as a model. Three applications cover most practice. As a residual estimator — the least-measurable term, almost always evapotranspiration, is obtained by difference from the measured ones; this is exactly Problem 7(i). As a continuous simulator — applied at a daily or hourly step with each store given a governing relation (an infiltration equation to move water from surface to soil, a recession to move it from groundwater to stream, a degree-day or energy-balance routine for snow), the balance is stepped forward to generate a continuous flow series; every conceptual rainfall–runoff model (HBV, HEC-HMS, HYDROTEL, GR4J) is this. As a design and operating tool — applied to a reservoir it gives yield, drawdown and refill; applied to a detention pond it gives the storage needed to hold a design event.

Worked example — monthly balance on a 250 km² catchment. Over a 30-day month the basin receives $P = 85$ mm, the outlet gauge averages $Q = 6.5\ \text{m}^3/\text{s}$, and soil-moisture probes indicate a storage decrease of 5 mm. Convert the discharge to a basin depth first:

$$R = \frac{Q\,\Delta t}{A} = \frac{6.5\ \text{m}^3/\text{s} \times (30 \times 86\,400)\ \text{s}}{250 \times 10^{6}\ \text{m}^2} = 0.06739\ \text{m} = 67.39\ \text{mm}$$

Then the residual, which lumps evapotranspiration with any net deep drainage, follows from the balance:

$$E + T + G = P - R - \Delta S = 85 - 67.39 - (-5) = \boxed{22.61\ \text{mm}}$$

That is a plausible cool-season ET for a Canadian catchment, and the check that it is plausible — positive, and consistent with the season and the vegetation — is the first validation any balance must pass. Note the unit discipline the calculation demands: discharge in m³/s must become a depth over the basin before it can be added to a rainfall in millimetres.

(ii) Three key steps in validating a watershed rainfall–runoff model (6 marks)

The conceptual diagram supplied with the question shows the model’s structure: climate input feeding a snow store and an evapotranspiration flux, a soil-moisture store beneath, and runoff generation followed by routing to the outlet. Validation must test that chain — not just the final hydrograph, because a model can produce the right outflow for the wrong internal reasons.

Step 1 — split-sample testing on an independent period. Divide the observed record in two. Calibrate the parameters on the first block, then run the model unchanged over the second block and compare simulated against observed discharge. Because the validation block was never used to fit anything, agreement there is evidence that the model has captured basin behaviour rather than memorised the calibration data. The test is made harder, and more informative, by a differential split sample — calibrating on wet years and validating on dry ones, or vice versa — which is the only way to probe whether the model will hold under the conditions it will actually be used to extrapolate to.

Step 2 — verify the internal state variables, not just the outflow. The diagram makes the reason explicit: outflow is the end of a chain, and compensating errors in the chain can still produce an acceptable hydrograph. Check the snow module against measured snow-water equivalent from snow-course or snow-pillow data, including the melt timing; check simulated evapotranspiration against a Penman–Monteith or eddy-covariance estimate and against seasonal norms; check soil moisture against probe data; check baseflow against the observed recession, and simulated recharge against water-table fluctuations. In a snowmelt-dominated Canadian basin this step is decisive, because a snow routine that melts a week early can be tuned to give the right annual volume while getting every freshet peak wrong.

Step 3 — apply objective goodness-of-fit criteria together with a bias and event check. Report quantitative measures — Nash–Sutcliffe efficiency, Kling–Gupta efficiency, coefficient of determination, RMSE — and, separately, the volume bias, since a model can score well on efficiency while systematically losing water. Then examine the events that matter: peak magnitude and timing for the largest floods, and the low-flow period, both of which a single aggregate statistic will hide. Finally, run a sensitivity and uncertainty analysis to identify which parameters actually control the result and to attach a confidence band to the output. Validation is not passed or failed by one number; it is a documented argument that the model is fit for the specific purpose it will be used for.

(iii) Muskingum routing and the continuity equation for reservoir and lake routing (6 marks)

Routing tracks a flood wave through a storage element — a reservoir, a lake, or a river reach — and answers what the outflow hydrograph looks like given the inflow. Every method starts from the storage form of continuity,

$$\frac{\mathrm{d}S}{\mathrm{d}t} = I(t) - O(t) \qquad\Longrightarrow\qquad \frac{S_2 - S_1}{\Delta t} = \frac{I_1 + I_2}{2} - \frac{O_1 + O_2}{2}$$

which is one equation in two unknowns ($S_2$ and $O_2$). A second, storage relation is needed to close it, and the choice of that relation is what distinguishes the methods.

Reservoir (level-pool) routing. In a reservoir or lake the water surface is essentially horizontal, so storage is a single-valued function of stage, and outflow through a spillway or outlet is also a single-valued function of stage. Storage therefore depends on outflow alone, $S = f(O)$, and continuity is solved by the storage-indication (modified Puls) method: rearranged as

$$\left(\frac{2S_2}{\Delta t} + O_2\right) = (I_1 + I_2) + \left(\frac{2S_1}{\Delta t} - O_1\right)$$

the right side is known at each step, and a pre-computed curve of $\left(2S/\Delta t + O\right)$ against $O$ — built from the elevation–storage and elevation–discharge relations — yields $O_2$ directly. This is the standard calculation for sizing a spillway to a design inflow flood or a stormwater detention pond, and it is the routing that Canadian Dam Association practice requires for the inflow design flood.

Muskingum (river/reach) routing. In a river reach the water surface is not horizontal during a flood: on the rising limb the inflow end stands higher than the outflow end, adding a wedge of storage above the level prism, and on the falling limb the wedge is negative. Muskingum represents this with

$$S = K\left[X I + (1 - X) O\right]$$

where $K$ is the reach travel time (the storage constant) and $X$ the weighting factor, $0 \le X \le 0.5$ — $X = 0$ recovers level-pool (reservoir) behaviour with maximum attenuation, and $X = 0.5$ gives pure translation with none. Natural channels typically fall between 0.1 and 0.3. Substituting into the discretised continuity equation gives the working form

$$O_2 = C_0 I_2 + C_1 I_1 + C_2 O_1$$

$$C_0 = \frac{\Delta t - 2KX}{2K(1-X) + \Delta t}, \quad C_1 = \frac{\Delta t + 2KX}{2K(1-X) + \Delta t}, \quad C_2 = \frac{2K(1-X) - \Delta t}{2K(1-X) + \Delta t}$$

with the essential check $C_0 + C_1 + C_2 = 1$, which guarantees mass conservation. For $K = 12$ h, $X = 0.20$ and $\Delta t = 6$ h the denominator is $2(12)(0.8) + 6 = 25.2$, giving

$$C_0 = \frac{6 - 4.8}{25.2} = 0.0476, \quad C_1 = \frac{6 + 4.8}{25.2} = 0.4286, \quad C_2 = \frac{19.2 - 6}{25.2} = 0.5238$$

and these sum to exactly 1.000. One routing step with $I_1 = 50$, $I_2 = 120$ and $O_1 = 45\ \text{m}^3/\text{s}$ gives

$$O_2 = 0.0476(120) + 0.4286(50) + 0.5238(45) = \boxed{50.7\ \text{m}^3/\text{s}}$$

The routing interval is not free: numerical stability and non-negative coefficients require $2KX \le \Delta t \le 2K(1-X)$, here $4.8 \le \Delta t \le 19.2$ h, and $\Delta t$ should also be short enough to define the rising limb (a common rule is $\Delta t \le t_p/5$). $K$ and $X$ are obtained from observed inflow–outflow pairs by plotting weighted discharge against storage and choosing the $X$ that collapses the loop to a straight line, or estimated from reach length and wave celerity where no record exists.