Continental Scale Heterogeneous Channel Flow Routing Strategy for Operational Forecasting Models

The benefits of operational forecasting models to the general public are numerous. They support water management decisions, provide the opportunity to mitigate the impacts of weather‐ and flood‐related disasters and potentially save lives and properties. Channel flow routing is a key component of these models and affects their ability to forecast flood depth, duration, and extent. Continental scale channel flow routing within the operational forecasting environments encounters a broad spectrum of hydraulic characteristics. Deploying computationally demanding approaches, such as the dynamic wave, should be limited in time and space to conditions where the inertia terms are significant (typically in low‐gradient environments and whenever backwater effects are prominent); otherwise, efficient and robust methods, e.g., Kinematic, Muskingum‐Cunge or diffusive waves should be the default. The heterogeneous routing approach presented here provides a framework to evaluate the balance between friction, inertia, and pressure and strategically triggers the appropriate wave approximation. The strategy recommended here is to activate the appropriate wave approximation based on the ambient hydraulic conditions, and smoothly transitions among these approximations. This strategy, if successfully implemented, would strike a balance among the performance metrics of operational forecasting models, namely, computational efficiency, accuracy, and minimization of computational instabilities.

length of the channel network of a given watershed. In turn, flood elevations and duration along the channel length affect the flood extent over the floodplains (Meselhe and Holly 1993;Meselhe et al. , 2004Sabur and Steffler 1996;Thielen et al. 2009;Pappenberger et al. 2011;Price et al. 2012;Burek et al. 2013).
Channel flow routing methods, to varying degrees, require information on the channel dimensions and longitudinal slope as well as the bed and bank roughness. Routing methods are based on numerical solutions of the well-known De St. Venant equations for mass and momentum as the primary governing equations for the flow of water in channels (Henderson 1966;Cunge et al. 1980). Applied in their full form, they are often labeled as the dynamic wave approximation, or "hydraulic routing." The dynamic wave approximation is computationally expensive relative to simpler approximation and requires a great deal of caution to avoid numerical instabilities (Cunge et al. 1980;Roberson et al. 1998). The computational expense and possible numerical instabilities are two concerning issues for operational models such as the National Water Model (https://water.noaa.gov/about/nwm).
If the inertial terms are omitted from the governing equations, the result is the diffusive wave approximation. If both the inertial terms as well as pressure difference terms are ignored, the result would be the kinematic wave approximation (Ponce 1989;Roberson et al. 1998). Together, the kinematic and diffusive wave routing methods are referred to as "hydrologic routing" or "bulk waves." These simplifications may result in substantial reduction in the computational load and enhancement of the stability and robustness. These are important consideration factors, especially for operational forecasting models. However, it should be noted that there is a tradeoff where the simplifications, while providing computational efficiency and robustness, may not be appropriate under certain hydraulic conditions. Channel flow routing is well documented in the literature (e.g., Weinmann and Laurenson 1979;Sabur and Steffler 1996;Munier et al. 2008). Channel flow routing can be performed using either dynamic, diffusive, or kinematic wave approximations. Another commonly used channel flow routing method is the Muskingum-Cunge method (Cunge 1969;Koussis 1976Koussis , 1980Ponce and Yevjevich 1978;Weinmann and Laurenson 1979;Smith 1980;Ponce 1983;U.S. Army Corps of Engineering 1991). Fundamentally, the Muskingum-Cunge method is a kinematic wave routing method along with some refinements that allow it to attain attributes of the diffusive wave equation (Smith 1980). As such, the Muskingum-Cunge method accounts for hydrograph convection and diffusion, which indicates that it captures some downstream hydrograph peak attenuation. It is computationally robust and has been thoroughly evaluated against field data and showed favorable agreement (U.S. Army Corps of Engineering 1991). It is more computationally efficient than the full diffusive wave approximation, and at times, it performs favorably and comparably against the diffusive, and even the dynamic, wave approximations. However, like other kinematic wave and diffusive wave routing methods, the Muskingum-Cunge method does not fully capture strong backwater effectsespecially in low-gradient (hydraulically mild slope) channels (U.S. Army Corps of Engineering 1991).
Of high relevance to the research presented here, it is vital to fully understand the range of applicability of various wave approximations.  and Ferrick (1985) provided a comprehensive overview and quantitative estimates of the range of validity of various wave approximations. It should be understood that a choice of a "universal" or a "homogenous" routing method would be adequate only within its range of applicability governed by the ambient hydraulic conditions. This is especially true for continental scale applications of operation and forecast models, where a broad spectrum of hydraulic conditions ranging from steep and flashy to low-gradient and attenuated responses are encountered. Thus, these applications may need to consider a "heterogeneous" routing method where the various wave approximations are deployed based on the ambient hydraulic conditions. Essentially, the robust and computationally efficient bulk wave routing approaches should be deployed as a default, whereas the full dynamic wave approximation should be used only when and where necessary. The approach to determine this "necessity" and an overall description of the heterogeneous routing strategy are outlined and discussed in the following section.

Governing Equations
The De St. Venant equations govern the conservation of flow volume and momentum through open channels, and are given by the following equation: where Q is the flow rate, A is the flow area, g is the gravitational acceleration, y is the water depth, S o is the bed slope, S f is the energy slope, x is the streamwise distance, and t is time. Equation (2) represents a balance between inertia, pressure, and friction. The appropriate wave approximation can be determined by that balance. It should be noted that friction is always important. Under certain ambient hydraulic conditions, inertia is significant, but overall rarely dominant (Ferrick 1985). When friction is dominant, Equation (2) can be approximated accurately by bulk waves (either the kinematic or diffusive). In the instances where inertia is significant, dynamic wave approximation would be required. Clearly, the ambient hydraulic conditions (the balance between friction and inertia) change temporally in response to rain events or the operation of water management structures. Spatially, the friction-inertia balance will also greatly vary among watersheds that are steep and those with low-gradient topography. The change in the spatiotemporal friction-inertia balance will be reflected through a change in the relative magnitude of the various terms in the momentum equation.
For the heterogeneous routing strategy to be effective in terms of selecting the appropriate wave approximation, triggers or flags must be identified. These triggers will be used to guide the selection of the wave approximation to be deployed. Clearly, calculating these triggers must be substantially less computationally expensive than calculating the full dynamic De St. Venant equations. Ferrick (1985Ferrick ( , 2005 suggested calculating dimensionless scaling parameters (DSP) as triggers or indicators of the appropriate wave approximation that needs to be deployed. The work presented here follows the same approach.

Heterogeneous Routing Strategy
The idea of the heterogeneous channel flow routing strategy revolves around determining the balance between friction and inertia to guide the selection of the appropriate wave approximation to be deployed in a given channel segment.
The work of Ferrick and Waldrop (1977), Ferrick (1979) and Ferrick et al. (1984) suggest that it is possible to identify the appropriate wave approximation by determining the wave celerity. When Equations (1 and 2) are differentiated with respect to x, then Equation (2) is differentiated with respect to t, the resulting equations can be combined into a single system equation (see Ferrick 1985 for details). This system equation can be written in dimensionless form to facilitate determining the friction-inertia balance as follows: where v* is v/v 0 , y* is y/y 0 , x* is x/Dx, and t* is t/Dt; whereas C r , F 0 , S, D I , F I , F c , and D are DSP defined in Table 1. where v 0 is the mean flow velocity; y 0 is the mean flow depth; Dx is the characteristic length scale of wavelength, Dt is the time required for the wave to travel distance Dx; Dx/Dt represents the measured wave celerity, and k is a water surface parameter (one for free surface; two for ice covered). The United States Geological Survey (USGS) devoted efforts to estimate river waves and the time of travel based on their gaging network (Jobson 2000). It has also been documented in the literature that bulk wave approximations are most suitable for long flood waves (Sriwongsitanon et al. 1998;Ferrick 2005;Humphries et al. 2014). Furthermore, the bulk wave approximations, as mentioned earlier, are prevalent when friction is dominant in the momentum equation. In such cases, the kinematic and diffusive approximations travel with a celerity governed by friction (related to the flow velocity) which is much slower than the dynamic wave celerity. The kinematic wave is further limited to channels with hydraulically steep slopes and where the channel reaches are divided into short segments. This is due to the fact that kinematic waves do not attenuate (while diffusive waves do).
It has also been discussed in the literature (Stoker 1957), and is broadly perceived, that dynamic waves are always present and therefore using the kinematic or the diffusive waves are less accurate. While this statement might be true, under certain hydraulic conditions the difference between the dynamic wave approximation and the bulk wave approximation might be negligible (especially for operational forecasting models). The momentum equation in all wave approximation cases remains nonlinear and is solved within computational models using some numerical approach. Numerical algorithms inherently include truncation errors that can lead to either diffusive or dispersive numerical error. Furthermore, when some Dimensionless diffusion coefficient: ratio of wave diffusion to wave advection terms (based on physical processes) are much smaller than others in the same set of equations, they become numerically "stiff" and are challenging to solve numerically. For example, if the inertia terms are extremely small under certain hydraulic conditions, attempting to solve the full dynamic equation for these conditions may require using higher order numerical methods to maintain stability. However, that could be too computational costly for operational forecasting models. For these applications especially at the continental scale, deploying the bulk wave approximations for hydraulic conditions where the inertia terms are small is advantageous as they provide computational speed while not sacrificing accuracy. Based on the analysis performed in Ferrick (1985), F c values larger than one (~10 or higher) indicates dominance of friction. In such cases, bulk wave approximations are an appropriate choice. In these cases, C r is also typically high (≥1). As the value of F c decreases (≤1) inertia becomes significant and the dynamic wave approximation may be necessary; C r also is in the order of, or lower than, 1. When F c or F I are much larger than 1 and concurrently the dimensionless diffusion coefficient, D, is less than 1, it indicates that advection dominates diffusion and hence, kinematic wave should be used. When F c or F I are much larger than 1, whereas D takes on a value of the order of 1, diffusive wave is most appropriate.
It should be emphasized that the transition among various wave approximations is not sudden; rather it is a continuous and gradual transition. Therefore, for a heterogeneous strategy to function efficiently, great care must be taken to avoid numerical instabilities or discrepancies in the solution at the transition zones.
The DSP, especially F c , F I , and D, can be used as indicators for the relative dominance of the momentum equations terms, they will be used as the foundation of the heterogeneous strategy to determine the selection of the wave approximation.
In the following section, we present a set of numerical experiments to explore how DSP can be used to determine the selection of the appropriate wave approximation.

ANALYSIS
The experiments discussed below were performed using two dynamic wave approximation channel flow routing models, namely MESH , and Preissmann Scheme coded in the HEC-RAS model (U.S. Army Corps of Engineering 2016). Being open source code, the MESH model allowed us to perform large number of permutations in a computationally efficient manner. Employing the publicly available HEC-RAS model allowed us to confirm that both dynamic wave solutions (Preissmann or MESH) produce consistent output.

Numerical Experiments Using the MESH Scheme
A set of numerical experiments were performed to examine the relationship between the DSP and the relative magnitude of various terms in the momentum equation (Equation 2). The experiments were performed using a dynamic wave approximation channel flow routing model, MESH . The geometrical attributes of the channel used in these experiments are summarized in Table 2. The channel dimensions, length, and bottom width, as well as spatial and temporal discretizations were kept constant for all experiments. We varied the following set of parameters: channel slope, the bed roughness coefficient, Manning's n (n), the volumetric flow rate (Q), and the downstream water depth. It should be noted that these experiments did not capture the full spectrum of all possible hydraulic conditions; rather we intended to demonstrate the relationship between the DSP and the routing wave approximations. We performed these experiments while varying one parameter at a time as shown in Table 3.
The flow hydrographs were generated using a skewed Gaussian distribution (see Figure 1), changing only the peak flow between experiments. It should be noted that the governing equations are numerically satisfied at the conclusion of individual time steps within the unsteady flow experiments, hence each time step could be treated as separate realization representing hydraulically valid instance of a flow profile.

Numerical Experiments Using the Preissmann Scheme (HEC-RAS)
The same simple channel described earlier was used again through this set of experiments using the Preissmann Scheme coded in the HEC-RAS model (U.S. Army Corps of Engineering 2016). We examined additional hydraulic conditions to supplement the experiments performed using the MESH scheme. The bed slope ranged from 10 À4 to 10 À3 and for some experiments we examined adverse slope of 5 9 10 4 . For some experiments we used normal depth as the downstream tailwater and for others we imposed a high water-level to induce backwater profiles. For all the experiments we used a triangular upstream flow hydrograph where the initial flow was 100 m 3 /s remaining constant for 30 min then linearly increasing to 600 m 3 /s within 2.5 h, then linearly decreasing to 100 m 3 /s within another 2.5 h and finally remaining constant at 100 m 3 /s until the end of the experiment.

Red River, Louisiana, USA Calculations
In addition to the idealized channels described above, a reach of the Red River in northwestern Louisiana was also analyzed (see Figure 2). The upstream location for this analysis was set to be the USGS river discharge observation gauge on the Red River near Hosston, Louisiana (USGS 2019a) and the downstream location was the USGS gauge on the Red River at Shreveport, Louisiana (USGS 2019b). The distance between these two USGS stations is 56 km, whereas the average bed slope is 7.7 9 10 À7 ; this value, along with the geometry of the channel cross sections at these two locations were determined from a HEC-RAS hydraulic model of the Red River which was downloaded from the FEMA Flood Risk Studies Engineering Library (FEMA 2019). This HEC-RAS model (FEMA 2011) was developed for a single reach of the Red River from the Louisiana-Arkansas border to a location 52 km downstream of Shreveport, Louisiana. Initially developed for a steady-state simulation of a so-called 100-year flood event, daily mean discharge data for the month of February 1966 recorded at the USGS gauge near Hosston, Louisiana (USGS 2019a) were used as upstream boundary condition for an unsteady flow simulation within the HEC-RAS model; a downstream boundary condition of normal flow depth was assumed. A single flood wave was modeled, which had a peak flow rate of 1,246 m 3 /s, which arrived on February 15, 1966. The peak discharge at the Shreveport gauge arrived one day later; only daily mean discharge data were available at the upstream and downstream locations. The HEC-RAS model  simulated the peak discharge arriving at Shreveport 24.5 h after it was at Hosston. This translates to a celerity of 0.64 m/s. Simulated water level and discharge were saved at a one-minute time step; HEC-RAS calculated stage-area and stage-conveyance curves were used to determine cross-sectional flow area and hydraulic radius at every one-minute data point at both the upstream and downstream cross sections. These hydraulic parameters and simulated hydraulics were then used to calculate the various  Range of applicability of the diffusive wave DSP to determine the appropriate wave approximation to be used. In addition to the Red River simulation described above with a normal depth as a downstream boundary condition, a second simulation was conducted with an artificially imposed backwater at the downstream location; a constant stage of 49.7 m was set at the furthest downstream cross section, which was 0.7 m above the simulated peak stage at the Shreveport gauge location when normal depth at the downstream boundary was assumed in the simulation. The details of the calculations for both of these cases are provided in the Supporting Information of this paper.

RESULTS AND DISCUSSIONS
The MESH dynamic wave model ) was used to perform the set of experiments described in the previous section. The momentum equation terms were calculated at each computational point and for each time step and were considered as individual realizations. A total of 67,724,800 realizations resulted from the experiments performed. Of that total, a 0.1% randomly selected sample was used to analyze the relative magnitude of the momentum equation terms as well as identify a relationship between selected DSP and the various river wave types. The sampling was done to reduce the computational effort needed and to facilitate plotting the results. We repeated the sampling process to ensure that it does not influence the trends or the conclusions. Figure 3 below shows the relative magnitude of the momentum equation terms. The inertia terms in Equation (2), namely local and convective accelerations, are referred to as Term III. The pressure gradient term in Equation (2) is Term II, whereas Term I refers to the friction and bed slopes. In essence, Equation (2) can be written as follows: Since friction is always important, we can normalize the equation by Term I as follows: III=I þ II=I þ 1 ¼ 0 or II=I ¼ À1ÀIII=I: It should be noted that when inertia is negligible, the ration II/I approaches À1. Table 4 provides a summary of the relative magnitude of the three terms calculated from the numerical tests in MESH and shown in Figure 3. It can be observed from Figure 3 and Table 4 that inertia terms are significant (more than 10% of the friction term) in approximately 24% of all the realizations. Furthermore, Table 4 clearly shows that where bulk waves are applicable, the diffusive wave approximation is almost always required (97%), leaving the kinematic wave with a narrow range of applicability. Now, we explore the relationship between F c and the applicability of the dynamic wave approximation. Per the analysis presented in Ferrick (1985), the dynamic wave approximation is needed when F c is in the order of 1. As F c approaches 10 (and higher) the bulk waves are appropriate. The pattern observed in Figures 4 and 5 is complex and requires careful discussion, but overall it shows that approximately 60%-80% of the time, bulk waves are applicable.
As discussed earlier, when bulk waves are applicable, a low value for the DSP D (~0.1) indicates a kinematic wave, where higher values indicate a diffusive wave. The variation of D vs. the relative magnitude of Terms I and II is shown in Figure 6. The ratio between Terms I and II show that they are predominantly in the same order of magnitude (and as illustrated in Figure 3 are of opposite sign) and hence the diffusive wave is almost always needed whenever bulk waves are applicable.
The statistics provided in Tables 4 and 5 summarize the trends presented in Figures 4-6. Both the tables and figures suggest that the dynamic wave approximation is needed for approximately 20%-30% of the calculated realizations. Where bulk waves are applicable, the diffusive wave is almost always needed leaving the kinematic wave with a narrow window of applicability.
The analysis presented in Figures 7-10 was performed using the HEC-RAS software (based on the Preissmann scheme). The analysis explores the relationship between the DSP, F c , F I , D, and C r and the ambient hydraulic conditions. The experiment presented in Figure 7 is for a relatively flat channel (bed slope of 10 À4 ) and where  The figure shows that F c and F I are both sufficiently high (~10) indicating a bulk wave applicability. The figure also shows a value of~1 for D indicating that a diffusive wave is needed for this application. This experiment indicates that even for relatively low-gradient channels if no backwater effects are experienced, a diffusive wave approximation is applicable. We repeated the same experiment presented in Figure 7 while imposing sufficiently high tailwater at the downstream end of the channel to induce backwater effects. This adjustment, shown in Figure 8, drastically changed the balance between inertia and friction. In this experiment, and as shown in Figure 8, F c and F I , representing friction, were in the order of~0.1 and 1, respectively, indicating the necessity of dynamic wave approximation.
We increased the bed slope from 10 À4 to 10 À3 while imposing a normal depth at the downstream end to ensure backwater effects are not present. For these conditions, as shown in Figure 9, both F c and F I are high (>10) indicating a bulk wave applicability, whereas the figure also shows a value of~0.1 or less for D indicating that a kinematic wave is applicable. Finally, and as expected, for the adverse slope channel, Figure 10, both F c and F I are low (clearly <10) indicating a dynamic wave applicability. Finally, the analysis we performed for the Red River is shown in Figure 11. The bed slope for this segment of the river is quite low, 7.7 9 10 À7 . A normal depth was imposed at the downstream end. As seen in Figure 11, and even for such a low-gradient channel, F c and F I are high (>10) indicating a bulk wave applicability. The scaling parameter D is in the order of 1 indicating the need for diffusive wave. This confirms the same behavior observed in Figure 8 that even for lowgradient channels, as long as backwater effects are not present, bulk waves are applicable. The effect of backwater on the scaling parameter F c for the Red River is shown in Figure 12. Compared to Figure 11 where normal depth was imposed at the downstream end, the value of F c dropped significantly when backwater effects were present ( Figure 12) indicating that the dynamic wave approximation might be needed.

CONCLUSIONS
Typically, the primary performance metrics for continental scale operational forecasting models are computational efficiency/speed; accuracy, and robustness (avoidance of computational failure). The channel flow routing component of these large-scale models affect their ability to capture flood depth, duration, and spatial extent. For continental scale application, the ambient hydraulic conditions; namely bed slope, roughness, presence of road-crossing, and other flow-impedance features, are strongly variable in time and space. To accommodate such variability, it is recommended to deploy a heterogeneous routing strategy.
The momentum equation articulates a balance among inertia, pressure gradient, and friction. This balance is governed by the ambient hydraulic conditions. The analysis we presented here illustrates that the appropriate river wave approximation is dictated by that balance and specifically by the relative significance of each of these processes. Physically, friction is always important and cannot be ignored. When friction is dominant, bulk waves are applicable and the utility of the dynamic wave approximation is not needed. Inertia, however, is rarely dominant, but is significant under certain circumstances.
The analysis we presented shows that the inertia terms, and consequently the dynamic wave approximation, are needed where backwater effects are prevailing. The larger the deviation from the uniform/ normal water depths the more significant the inertia terms become necessitating the deployment of the dynamic wave approximation. If no backwater effects are present, and even in low-gradient channels, the bulk waves are applicable. In these instances, the diffusive wave is almost always needed leaving the kinematic wave with a limited range of applicability.
Using the DSP, F c , F I , and D, is an effective approach to identify and trigger the appropriate wave approximation. When F c (or F I ) is sufficiently large (much larger than 1) the bulk waves (diffusive or kinematic) are applicable, whereas for lower values (in the order of 1) the dynamic wave approximation is needed. When bulk waves are applicable, the . Time evolution of three scaling parameters; F c , F I , and C r for the Red River for the backwater downstream boundary condition case. The bed slope for this segment is 7.7 9 10 À7 . kinematic wave approximation is appropriate when D is much smaller than one; otherwise the diffusive wave approximation would be needed. Future research would focus on applying the three wave (kinematic, diffusive, and dynamic) approximations, independently, to cases where data are available. We will quantify the error resulting from each wave approximation to carefully delineate the values of the scaling parameters that should be implemented in a heterogeneous routing approach. We will also quantify the computational cost of each of the three wave approximations to determine the tradeoffs while delineating the transitions among the wave approximations.

SUPPORTING INFORMATION
Additional supporting information may be found online under the Supporting Information tab for this article: Excel file showing the details of the calculations and experiments discussed in the manuscript.