Modified Puls Model

Basic Concepts and Equations

The Modified Puls routing method, also known as storage routing or level-pool routing, is based upon a finite difference approximation of the continuity equation, coupled with an empirical representation of the momentum equation (Chow, 1964; Henderson, 1966). For the Modified Puls model, the continuity equation is written as

This simplification assumes that the lateral inflow is insignificant, and it allows width to change with respect to location. Rearranging this equation and incorporating a finite-difference approximation for the partial derivatives yields:

where $//<![CDATA[ \begin{array}{l}\overline{I_{t}}\end{array} //]]>$ = average upstream flow (inflow to reach) during a period $//<![CDATA[ \begin{array}{l}\Delta t\end{array} //]]>$; $//<![CDATA[ \begin{array}{l}\overline{O_t}\end{array} //]]>$ = average downstream flow (outflow from reach) during the same period; and $//<![CDATA[ \begin{array}{l}\Delta S_t\end{array} //]]>$ = change in storage in the reach during the period. Using a simple backward differencing scheme and rearranging the result to isolate the unknown values yields:

in which $//<![CDATA[ \begin{array}{l}I_{t-1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}I_t\end{array} //]]>$ = inflow hydrograph ordinates at times t-1 and t, respectively; $//<![CDATA[ \begin{array}{l}O_{t-1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}O_t\end{array} //]]>$ = outflow hydrograph ordinates at times t-1 and t, respectively; and $//<![CDATA[ \begin{array}{l}S_{t-1}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}S_t\end{array} //]]>$ = storage in reach at times t-1 and t, respectively. At time t, all terms on the right-hand side of this equation are known, and terms on the left-hand side are to be found. Thus, the equation has two unknowns at time t: $//<![CDATA[ \begin{array}{l}S_t\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}O_t\end{array} //]]>$. A functional relationship between storage and outflow is required to solve Equation 3. Once that function is established, it is substituted into Equation 3, reducing the equation to a nonlinear equation with a single unknown, $//<![CDATA[ \begin{array}{l}O_t\end{array} //]]>$. This equation is solved recursively by the program, using a trial-and-error procedure. Note that at the first time t, the outflow at time t-1 must be specified to permit recursive solution of the equation; this outflow is the initial outflow condition for the storage routing model.

If the storage vs. discharge relationships are carefully constructed using a hydraulic model that includes bridges and/or other hydraulic structures, this method can simulate backwater effects and the impacts of hydraulic structures so long as the effects/impacts are fully contained within the reach.

Required Parameters

The parameters that are required to utilize this method within HEC-HMS are the initial condition, a storage vs. discharge relationship, and the number of subreaches.  An optional elevation vs. discharge function can be selected in addition to an optional invert.

Two options for specifying the initial condition are included: outflow equals inflow and specified discharge [ft³/sec or m³/sec].  The first option assumes that the initial outflow is the same as the initial inflow to the reach from the upstream elements which is equivalent to the assumption of a steady-state initial condition.  The second option is most appropriate when there is observed streamflow data at the end of the reach.  In either case, the initial storage will be computed from the first inflow to the reach and corresponding storage vs. discharge function.

The storage vs. discharge relationship (which is a paired data object) must define the amount of outflow for a specific amount of storage in the reach.  Entered storage values should cover the entire range of expected storages that will be encountered during a simulation; the first storage value should be set to zero.  Typically, hydraulic routing simulations that include bridges and/or hydraulic structures are used to generate the storage vs. discharge relationships (Hydrologic Engineering Center, 2023).

The specified number of subreaches affects attenuation.  One subreach results in the maximum amount of attenuation and increasing the number of subreaches approaches zero attenuation.  For an idealized channel, the travel time through a subreach should be approximately equal to the simulation time step.  An initial estimate of this parameter can be obtained by dividing the actual reach length by the product of the wave celerity and the simulation time step.  For natural channels that vary in cross-section dimension, slope, and storage, the number of subreaches can be treated as a calibration parameter.  The number of subreaches may be used to introduce numerical attenuation which can be used to better represent the movement of flood waves through the natural system.

The optional elevation vs. discharge function should represent the depth of water for any given outflow from the reach.  When an elevation vs. discharge function is specified, the optional invert elevation should also be specified.

Defining the Storage-Outflow Relationship

The storage-outflow relationship required for the Modified Puls routing model can be determined using several techniques.

Using Hydraulic Model Outputs

Water-surface profiles can be computed with a hydraulic model for a range of discharges.  Hydraulic modeling applications like HEC-RAS (USACE, 2023) include automated tools that can be used to define a relationship of storage to flow between two channel cross sections using computed water-surface profiles.

The following figure illustrates a set of water-surface profiles between cross section A and cross section B of a channel. These profiles were computed for a set of steady flows, $//<![CDATA[ \begin{array}{l}Q_1\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}Q_2\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}Q_3\end{array} //]]>$, and $//<![CDATA[ \begin{array}{l}Q_4\end{array} //]]>$. For each profile, the volume of water in the reach, $//<![CDATA[ \begin{array}{l}S_i\end{array} //]]>$, can be computed, using solid geometry principles. In the simplest case, if the profile is approximately planar, the volume can be computed by multiplying the average cross-section area bounded by the water surface by the reach length. Otherwise, another numerical integration method can be used. If each computed volume is associated with the steady flow with which the profile is computed, the result is a set of points on the required storage-outflow relationship. This procedure can be used with existing or with proposed channel configurations. For example, to evaluate the impact of a proposed channel project, the channel cross sections can be modified, water surface profiles recalculated, and a revised storage-outflow relationship developed.

Using Historical Observations

Storage-outflow relationships can also be determined using historical observations of flow and stage. Observed water surface profiles, obtained from high water marks, can be used to define the required storage-outflow relationships, in much the same manner that computed water-surface profiles are used. Each observed discharge-elevation pair provides information for establishing a point of the relationship.  Sufficient stage data over a range of floods is required to establish the storage-outflow relationship in this manner. If only a limited set of observations is available, these values may be better suited to calibrate a water-surface profile-model for the channel reach of interest. Then the calibrated model can be used to establish the storage-outflow relationship as described above.

Trial and Error

Finally, storage-outflow relationships can be calibrated using observed inflow and outflow hydrographs for the reach of interest. Observed inflow and outflow hydrographs can be used to compute channel storage by an inverse process of flood routing. When both inflow and outflow are known, the change in storage can be computed using Equation 3. Then, the storage-outflow function can be developed empirically. Note that tributary inflow, if any, must also be accounted for in this calculation.  Inflow and outflow hydrographs also can be used to find the storage-outflow function by trial-and-error. In that case, a candidate function is defined and used to route the inflow hydrograph. The computed outflow hydrograph is compared with the observed hydrograph. If the match is not adequate, the function is adjusted, and the process is repeated.

Estimating Other Model Parameters

The Kinematic Wave Transform Model section of this manual describes how an accurate solution of the finite difference form of the kinematic wave model requires careful selection of $//<![CDATA[ \begin{array}{l}\Delta x\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\Delta t\end{array} //]]>$; this is also true for solution of the storage-routing model equations. For the kinematic wave model, an accurate solution can be found with a stable algorithm when $//<![CDATA[ \begin{array}{l}\Delta x/\Delta t \approx c\end{array} //]]>$ , where c = average wave speed over a distance increment $//<![CDATA[ \begin{array}{l}\Delta x\end{array} //]]>$. This rule also applies with storage routing. As implemented in the program, $//<![CDATA[ \begin{array}{l}\Delta x\end{array} //]]>$ for the finite difference approximation of $//<![CDATA[ \begin{array}{l}\partial Q/\partial x\end{array} //]]>$ is implicitly equal to the channel reach length, L, divided by an integer number of steps. The goal is to select the number of steps so that the travel time through the reach is approximately equal the time step $//<![CDATA[ \begin{array}{l}\Delta t\end{array} //]]>$. This is given approximately by:

The number of steps affects the computed attenuation of the hydrograph. As the number of routing steps increases, the amount of attenuation decreases. The maximum attenuation corresponds to one step; this is used commonly for routing though ponds, lakes, wide, flat floodplains, and channels in which the flow is heavily controlled by downstream conditions. Strelkoff (1980) suggests that for locally-controlled flow, typical of steeper channels:

where y_o = normal depth associated with baseflow in the channel. Engineer Manual 1110-2-1417 Flood-Runoff Analysis (U.S. Army Corps of Engineers, 1994) indicates that this parameter, however, is best determined by calibration, using observed inflow and outflow hydrographs.