Following is a sample equation for Outflow:

Following is the equation for flow fraction:

This method can also solve for inflow, given outflow.

This method executes when the Inflow is known for the current timestep. The Inflow value is then used to calculate the Variable Lag Time. Outflow values can then be solved for at the timesteps corresponding to the Variable Lag Time. Following is a sample calculation:

In these equations, the integer value of a number means that the number is truncated after the decimal point with no rounding.

Since the lag time can vary with every timestep, it possible that more the one Inflow can contribute to a particular Outflow. Return Flow and Total Gain Loss can then be added to the calculated Outflow value.

There must be the same number of values in the Lag Coeff table as the value given in Num of Coeff in order to successfully dispatch at any timestep. If any needed value is invalid, an error will occur, and the run will stop.

The general equation for this method is as follows:

The general calculations for this method are very similar to the Impulse Response routing method.

This Routing method is executed from the Solve given Inflow dispatch method; that is, it only solves for the Outflow given the Inflows, not vice versa.

In the following, Inflow refers to the sum of the Inflow and Local Inflow slot. The variable ncoeff is the number of rows in the Variable Lag Coefficients slot. The outflows will be calculated using the following algorithm:

Each C coefficient is selected from the appropriate column of the Variable Lag Coefficients slot based on the value of the Inflow with which it is multiplied. For a particular evaluation of the equation, the sum of the coefficients does not necessarily equal 1.0. But, mass is preserved over time. Using the above sample table, if the Inflow at the current timestep t, is 12,000cfs and the Inflow at t-1 is 8,000cfs, C0 is 0.30 and C1 is 0.35 using the above table.

If the method is looking for an inflow for a timestep that is actually past the current timestep and the value is not valid, this inflow will be assumed to be zero. If the inflow is not valid for a previous timestep, the calculation will exit the method without calculating an outflow, post a warning message, and move on to the next timestep.

Following is the general equation for this method:

Following is the general equation for this method:

You can further discretize the reach into n sub-reaches for better control over attenuation.

First, the method checks the following:

The solution is a nonlinear numerical approximation to the following continuity and momentum equations.

For this numerical solution, accuracy increases as the Courant number approaches unity and decreases as the Courant number diverges from unity. The Courant number for each numerical grid section for the current timestep can be observed in the Numerical Parameters Output under the column labeled Courant. The Courant number is calculated as follows:

This hydraulic routing method also requires the user to select from another category of user selectable methods, Depth to Flow. The inputs for these methods are described under that user method category.

In the following method, a superscript refers to the timestep, and a subscript refers to the spatial location (element number) within the reach, upstream to downstream.

The solution is a nonlinear numerical approximation to the following continuity and momentum equations.

For this numerical solution, the total mass balance error will decrease as the Delta X Computational Element Length (x) decreases (number of segments increases). Decreasing x will reduce the dispersion in the numerical solution, resulting in higher peaks and lower valleys in the Reach Outflow. Decreasing x will also increase the Courant number. The Courant number for each numerical grid section for the current timestep can be observed in the Numerical Parameters Output under the column labeled Courant. The Courant number is calculated as follows:

This hydraulic routing method also requires the user to select from another category of user selectable methods, Depth to Flow. The inputs for these methods are described under that user method category.

In the method description below, a superscript refers to the timestep, and a subscript refers to the spatial location (element number) within the reach, upstream to downstream.

At the start of the run this method sets all columns of the Distributed Flow Output at the initial timestep (Start Timestep - 1) equal to the initial Reach Inflow (that is, the flow is the same at every element). It then calls the selected Depth to Flow method to calculate all flow parameters for each element at the initial timestep.

For each timestep within the run, it first sets the flow for the first element, as follows:

If Diversion or Local Inflow are not used, they are set to zero.

Then the selected Depth to Flow method is called calculate all flow parameter values for the first element using . The values are set in the first columns of the Distributed Output table series slots.

The method then loops over all remaining elements in the reach, upstream to downstream, to carry out the finite difference approximation in the following steps.

(This equation is a rearrangement of equation 9.6.6 in Chow et. al., 1988². Celerity is substituted using equation 9.3.11.)

After calculating flow at the final element, all gains, losses and Return Flow are added to the final element flow to give the total Reach Outflow.

The standard form of the Muskingum equation is as follows:

When applied to a spatially distributed grid, the flow in grid i is as follows:

C and D are calculated as follows:

Q is the grid’s flow, So is the Energy Slope, c is the wave celerity, and Tw is the top width of the channel calculated for the given flow Q. The wave celerity c and Tw are found using the selected Depth to Flow method for the flow, Q, at the appropriate timestep (see Depth to Flow for details). The incremental timestep is specified by the user while the spatial distribution is calculated to keep the Courant number close to 1 for the reference discharge.

The spatial distribution is automatically calculated by RiverWare. The incremental spatial step is determined such that the Courant number, C, will be close to one to reduce the effects of numerical dispersion. Since the discharge will vary for a simulation, the Courant number will also vary. To pick a value for the incremental spatial step that minimizes the effects of numerical dispersion, the user inputs maximum and minimum discharges expected for a simulation in the Extreme Values slot. The average of these two discharges is used as the reference discharge. The wave celerity, c, is then computed with this reference discharge using the selected Depth to flow; see Depth to Flow for details. The spatial step size, is computed with the input and the Courant number, C set to be 1.0.

The scheme is unconditionally stable for 0.0 <= X <= 0.5 or with D (cell Reynolds number) values less than unity. The cell Reynolds number, D, is defined as the ratio of hydraulic diffusivity to grid diffusivity, as follows:

The value of D for each numerical grid section for the current timestep can be observed in the Distributed MuskingumCunge Output table under the column labeled D.

The value of C + D for each numerical grid section for the current timestep can be observed in the Distributed MuskingumCunge Output table under the column labeled C + D.

The method proceeds as follows.

The standard form of the Muskingum equation is as follows:

When applies to a spatially distributed grid, the flow in grid i is as follows:

C and D are calculated as follows:

Q is the grid’s flow, So is the Energy Slope, c is the wave celerity, and Tw is the top width of the channel calculated for the given flow Q. The wave celerity c and Tw are found using the selected Depth to Flow method for the flow, Q, at the appropriate timestep (see Depth to Flow). The incremental timestep is specified by the user while the spatial distribution is calculated to keep the Courant number close to 1 for the reference discharge.

In the variable parameter, Muskingum-Cunge routing method, the cell Reynolds number, D, and the Courant number, C change at each timestep for each cell. Changing C and D values in the middle of a simulation results in a volume conservation error and storage inconsistency at steady state as shown by Todini. See Muskingum-Cunge for details.

To mitigate this loss of mass, Todini rederived the equations (not shown) and calculated the coefficients as follows:

The Muskingum-Cunge Improved method is an implementation of this scheme. Although Todini proposed relationships to define the channel parameters for common channel geometry, this method uses an iterative 3 and 4 point scheme and the selected Depth to Flow method to determine these same parameters (area, top width, depth, velocity, and celerity).

In addition, the reach inflow is used as the inflow to each incremental timestep. The Outflow is an average of the outflow over each incremental timestep. Full details of this method follow the slot descriptions.

The incremental timestep is related to the spatial distribution as follows:

Since the discharge varies for a simulation, the Courant number will also vary. To pick a value for the incremental spatial step that minimizes the effects of numerical dispersion, the user inputs minimum and maximum discharges expected for a simulation in the Extreme Flow Values slot. The average of these two discharges is used as the reference discharge when computing the wave celerity, c, using the selected Depth to Flow method. Then, the incremental spatial step is determined with the Courant number, C, set equal to 1.0 to reduce the effects of numerical dispersion.

The method proceeds as follows.

This method was implemented to achieve better mass conservation but is not guaranteed to fully conserve mass. Because the user selects the Extreme Flow values and Routing Timestep, there are circumstances where flow is not fully conserved. This is especially true when there are abrupt flow changes and flows that are zero. Testing has shown that mass is fully conserved in well behaved problems but not necessarily in real-world problems. Regardless, the mass conservation is much better than the original routing method.

The solution is a nonlinear numerical approximation to the following continuity and momentum equations.

For this numerical solution, accuracy increases as the Courant number approaches unity and decreases as the Courant number diverges from unity. The Courant number for each numerical grid section for the current timestep can be observed in the Numerical Parameters Output under the column labeled Courant. The Courant number is calculated as follows:

This hydraulic routing method also requires the user to select from another category of user selectable methods: Depth to Flow. The inputs for these methods are described under that user method category.

The algorithm uses is a predictor corrector scheme that is unstable for a Courant number greater than unity. The advantage of this scheme is that it is second order accurate and can minimize numerical diffusion. The general equations for this method are as follows.

The algorithm proceeds as follows.

The algorithm proceeds as follows.

First, the method checks the following.