Open-Channel Computation with R

The rivr package provides a computational toolset for simulating steady and unsteady one-dimensional ﬂows in open channels. It is designed primarily for use by instructors of undergraduate-and graduate-level open-channel hydrodynamics courses in such diverse ﬁelds as river engineering, physical geography and geophysics. The governing equations used to describe open-channel ﬂows are brieﬂy presented, followed by example applications. These include the computation of gradually-varied ﬂows and two examples of unsteady ﬂows in channels—namely, the tracking of the evolution of a ﬂood wave in a channel and the prediction of extreme variation in the water-surface proﬁle that results when a sluice gate is abruptly closed. Model results for the unsteady ﬂow examples are validated against standard benchmarks. The article concludes with a discussion of potential modiﬁcations and extensions to the package.


Introduction
Hydraulic engineers routinely rely on computational models to simulate the flow and conveyance of water for planning and designing channels and in-stream structures, managing reservoirs, and conducting flood damage and risk analyses. The simulation of fluid flow is a difficult and computationallyintensive endeavor, and simplifications are often made in order to reduce the task to one that is feasible with regard to time management and data availability. One of the most frequently-made assumptions is that of one-dimensional flow. This greatly simplifies the computation of a wide range of practicallyrelevant flows such as e.g. the prediction of the movement of a flood wave in a waterway, or the consequences of a sudden opening or closure of a control gate on the water level in a navigation channel. In making this assumption, the flow velocity-which in reality can vary in all three coordinate directions-is both depth-and laterally-averaged to produce a bulk or mean channel velocity that is a function of only time and distance along the channel: u = y=H y=0 z=S z=0 u(x, y, z, t) dy dz dy dz where u is velocity, y is the vertical dimension, z is the lateral (cross-channel) dimension, t is time and x is the along-channel direction. The integral bounds H and S are the flow depth and span, respectively. Treatment of bed friction is also frequently simplified in practical applications by relating the bed shear stresses to the "friction slope", which is derived from the Energy Grade Line (EGL) equation (Leopold, 1953) and expressed using the semi-empirical Manning's equation where S f is the friction slope, n is the empirical Manning's coefficient and R is the hydraulic radius. Manning's equation and the EGL equation are widely used in the discipline to predict steady and gradually-varied flow conditions. The prediction of unsteady flows is somewhat more complex; in addition to the fact that the flow area is different from the area of the conveying channel, the governing equations are both non-linear and coupled. This complicates the task of obtaining solutions that are both numerically stable and accurate without the need for excessive computational resources.
A wide variety of software suites are available for computing steady and unsteady flows; some of the most well-known packages are listed in Table 1. These packages are typically accessed via a graphical user interface and often have specialized input and output data formatting requirements. These software suites are capable of modeling complex networks of irregular channels and are therefore of great utility to engineers and practitioners working on large-scale projects, but they also pose high barriers to entry in terms of cost, usability and computer hardware requirements as well as in terms of ancillary data needs. These software suites are so complex that many developers provide training courses for utilizing the software, with the assumption that users already possess an education or background in river mechanics and hydraulics.
The rivr package (Koohafkan, 2015) was designed to provide rapid solutions to steady and unsteady one-dimensional flows in open channels. The primary purpose of the package is to provide educators with an accessible open-source toolset that can be used to supplement lectures with interactive examples, and to facilitate the design of assignments that encourage students to explore open-channel flow concepts rather than focusing on the actual implementation of solutions. Project for Statistical Computing provides an ideal platform for providing rivr (and educational tools in general) to students and instructors. R is both cross-platform and free and open-source software, meaning students can run the software on personal computers and educational institutions do not need enter expensive licensing agreements to provide it on institution-owned computers. Moreover, R provides a large and comprehensive library of functions for data visualization, analysis and export both through its core functionality and through a diversity of packages developed and supported by an active global community of over 2 million users, developers and scientists. R is often referred to as a "scripting language", meaning that code blocks or snippets can be run either in an automated fashion or piecewise by users; scripting languages are ideal for teaching because they support interactive programming and exploration of data. Furthermore, support for dynamic document generation with R via packages such as knitr (Xie, 2014) provides an excellent mechanism for instructors to develop teaching materials that place concept and implementation side by side, as well as a promising format for students to submit assignments and reports. Finally, R provides the potential for instructors with virtually no web development proficiency to create interactive web-based tutorials and demonstrations using the shiny web application framework (Chang et al., 2015).
The flow algorithms provided by rivr are implemented in C++ and integrated into the package using Rcpp (Eddelbuettel and Francois, 2011) to leverage the speed advantage of compiled code. All algorithms were developed from the ground up and are based on well-established, peer-reviewed contributions to the field. Functions are provided for computing normal and critical depth, channel geometry and conveyance, and gradually-varied flow profiles. The package also provides an advanced function which implements two canonical formulations of unsteady flow-the Kinematic Wave Model and the Dynamic Wave Model-which are routinely taught in graduate-level courses in hydrology and open-channel hydraulics. The numerical schemes implemented in rivr are similarly well-established and reflect varying degrees of sophistication that are common topics in numerical methods courses, including implementations of iterative solutions to implicit equations; a solution scheme for nonlinear ordinary differential equations; and numerical methods for solving sets of partial differential equations ranging from one of the simplest formulations of an explicit finite-difference scheme to an advanced two-step predictor-corrector scheme that illustrates the trade-offs between performance and complexity of implementation. Detailed explanation of the theoretical foundations and numerical schemes is provided in the form of a technical vignette, and the source code is carefully organized and commented to facilitate modifications and extensions to the package such as new numerical schemes and alternative formulations for steady and unsteady flow modeling.
This article provides use examples for computing both gradually-varied and unsteady flows with rivr. In the first section, solutions to gradually-varied flow or "backwater curve" problems using the standard-step method are discussed. In the second section, solutions to unsteady flow problems using different theoretical formulations are described, followed by a comparison of two numerical solutions to an unsteady flow problem with discontinuous boundary conditions. The article concludes with a discussion of potential opportunities for extensions and contributions to the package.

Calculation of back-water curves
The water-surface profile resulting from non-uniform steady flow conditions is generally referred to as a gradually-varied flow (GVF) profile or "backwater curve". Backwater curves are commonly used for floodplain delineation and assessing upstream effects of bridges, weirs and other channel structures. The GVF concept is based on a simplified form of the St. Venant equations (discussed in the following section) and is expressed as where dy/dx is the water-surface slope, S 0 is the channel slope and Fr is the Froude Number. Eq. (3) describes steady-state conditions, i.e. the channel flow rate Q = uA is constant. When dy/dx = 0, the flow depth is uniform and S 0 = S f ; the flow depth under these conditions is referred to as the The R Journal Vol. 7/2, December 2015 ISSN 2073-4859 "normal depth" and can be calculated using Manning's equation. The relative positioning of the normal depth and the "critical depth"-the depth at which the flow energy is minimized-determine the flow behavior. GVF profiles are classified based on the flow regime, the channel slope and the nature of the divergence from the steady-flow water-surface profile at a specific location or "control section" (Doubt, 1971). An M1 profile refers to sub-critical flows where the water depth at the control section is greater than the normal depth. An M2 profile refers to a water-surface profile where the water depth at the control section is below the normal depth, but above the critical depth. Other classes include "S" (steep channels where the normal depth is shallower than the critical depth), "C" (flow regimes where the normal depth and critical depth coincide), "H" (flat channels) and "A" (channels of "adverse" slope, i.e. sloped against the direction of flow).
The rivr package computes the GVF profile using the standard-step method (Chaudhry, 2007). A Newton-Raphson scheme is used where iterative methods are required to find the roots of an equation. The following example illustrates how to use rivr to calculate the variation of the water-surface profile with distance, and its return to the normal-depth level far upstream of a perturbation in flow depth. Both an M1 and an M2 profile are considered for a rectangular channel with a width of 100 feet, a slope of 0.001, a flow of 250 cubic feet per second, and a Manning's roughness coefficient of 0.045. The water depth y at the control section is specified as 1 foot above the normal depth y n for the M1 profile (2.71 feet) and 1.1 times the critical depth y c for the M2 profile (0.64 feet). The elevation and x-coordinate of the control section are both specified as zero. Example code for computing the water-surface profile using the function compute_profile is shown below. # calculate normal depth for channel, initial guess of 2 feet yn <-normal_depth(slope, mannings, flow, 2, Cm, width, sideslope) # calculate critical depth, initial guess of 1 ft yc <-critical_depth(flow, 1, g, width, sideslope) # compute the M1 profile 3000 ft upstream with a step size of 50 ft m1 <-compute_profile(slope, mannings, flow, yn + 1, Cm, g, width, sideslope, z0 = 0, x0 = 0, stepdist = 50, totaldist = 3000) # compute the M2 profile under the same assumptions m2 <-compute_profile(slope, mannings, flow, 1.1 * yc, Cm, g, width, sideslope, z0 = 0, x0 = 0, stepdist = 50, totaldist = 3000) The effect of model resolution on results is shown in Figure 1. For the M1 profile, the solution is quite stable even for very coarse (> 100 feet) resolutions. In contrast, solutions to the M2 profile are quite sensitive to model resolution owing to the steep gradient in dy/dx where y → y c . However, the solution to the M2 profile is not sensitive to resolutions finer than 10 feet. A step size of 10 feet is used in subsequent analyses to ensure smooth calculations of the water-surface slope for both profiles.
Manning's coefficient is an empirical quantity that cannot be directly measured, and engineers often rely on tabulated values to inform their selection of a roughness value based on channel bank material and conditions such as vegetation and compaction. It is generally necessary to investigate the sensitivity of the solution to the choice of Manning's roughness coefficient, or to calibrate a GVF model based on channel depth observations. The next example considers the GVF profile of the control section specified above for a range of Manning's roughness values spanning 0.0225 to 0.0675. Figure 2 shows that the normal depth (apparent as the section where the water-surface slope is linear) changes depending on the value of the roughness coefficient; this comes as no surprise, as the concept of normal depth is itself derived from Manning's equation. The specification of a fixed control section depth therefore confounds analysis of the relationship between Manning's coefficient and the nature of the GVF profile.
In order to control for the dependence of the normal depth on Manning's coefficient, the analysis is repeated using a control section depth of 1.25y n for the M1 profile and 0.75y n for the M2 profile, where y n is the normal depth associated with each roughness value. Note that the critical depth is by definition independent of Manning's coefficient. Figure 3 shows the water-surface elevation profile as a percent difference from the normal depth for each roughness value; it is clear from the bed roughness affects the distance that a depth disturbance propagates upstream under sub-critical flow conditions. Note that the M1 and M2 profiles are not symmetric; that is, the point at which the water depth returns to the normal depth from a 25% increase in water depth at the control section is further upstream compared to that for a 25% decrease in depth. This is because the velocity term in the energy equation is non-linear (u 2 ) and therefore the energy gradient is steeper for the M2 profile.

Computation of unsteady flows
Open-channel unsteady flows are typically described using the St. Venant equations, also known as the shallow water equations. The St. Venant equations are derived from the Navier-Stokes equations based on a number of simplifying assumptions, including uniform density and hydrostatic pressure conditions. The one-dimensional St. Venant equations are written as ∂A ∂t ∂Q ∂t where Eq. (4) is the continuity, or mass conservation, equation (assuming no infiltration or precipitation) and Eq. (5) is the momentum equation. The terms Q, u, A and S 0 and S f are as defined previously and g is gravitational acceleration. The term y results from averaging the pressure component and refers to the distance from the free surface to the centroid of the cross section. The rivr package provides multiple methods for solving these equations to compute unsteady flows through a prismatic channel. These methods are discussed in the following sections.

Simulations with the Kinematic Wave Model
The Kinematic Wave Model (KWM) provides a simplified solution to the St. Venant equations based on a truncated momentum equation that ignores inertial terms (Lighthill and Whitham, 1955). Eq. (5) is instead expressed through the relation where P is the wetted perimeter and S 0 is substituted for S f in Manning's equation.
The KWM engine (and other engines, as we show later) provided by rivr can be called via the function route_wave. The behavior of the KWM is demonstrated by using it to calculate the evolution of a flood wave down a 100-foot wide rectangular channel that is 150,000 feet long, with a slope of 0.001 and a Manning's roughness of 0.045. The boundary.condition argument is a vector specifying the flow at the upstream boundary for an arbitrary period of time in intervals equal to the timestep argument, in seconds. The monitor.nodes and monitor.times arguments are used to specify the function outputs. The hydrograph (change in flow over time) at the upstream boundary, shown in Figure 4, is modeled as a symmetrical flood wave and is conveniently expressed as where t is in seconds. # time steps to monitor initflow <-250 # initial flow condition (ft^3/s) times <-seq(0, 76000, by = dt) boundary <-ifelse(times < 9000, # upstream hyrograph (ft^3/s) 250 + (750/pi) * (1 -cos(pi * times/(4500))), 250) route_wave(slope, mannings, Cm, g, width, sideslope, initflow, boundary, The R Journal Vol. 7/2, December 2015 ISSN 2073-4859 timestep = dt, spacestep = dx, numnodes = numnodes, monitor.nodes = monpoints, monitor.times = montimes, engine = "Kinematic") While the KWM is easy to use and implement, it ignores wave damping and attenuation resulting in unrealistic predictions of the shape of the flood wave as it progresses downstream. The truncation of Eq. (5) results in what is known as "kinematic shock", i.e. the shape of the flood wave becomes increasingly asymmetrical as it progresses downstream. This effect is apparent in Figure 5, which shows the wave front approaching a vertical line as the flood wave travels along the channel.
The rivr package implements the KWM with an explicit, first-order accurate backwards-difference scheme. While such schemes are simple to conceptualize and implement, they are sensitive to numerical error and can be unstable. A necessary condition for stability is that the Courant Number C is less than unity, i.e.
where u is defined here as the initial velocity at the upstream boundary and ∆x and ∆t are the spatial and temporal resolution of the model. Numerical error is propagated through successive calculations and results in a systematic decrease in the accuracy of the solution as the flood wave progresses downstream, apparent in Figure 5 as a decreasing peak flow magnitude. The sensitivity of the solution generated by the KWM engine to model resolution is assessed by varying ∆x and ∆t while maintaining a Courant number of 0.7 to ensure numerical stability. The effect of numerical error on model accuracy is inferred by evaluating the shape of the hydrograph at a cross-section 50,000 feet downstream. Figure 6 shows the shape of the flood wave at said cross-section for selected spatial resolutions. As resolution decreases, the peak flow downstream decreases due to the effects of numerical error. In addition, the severity of the kinematic shock decreases with decreasing resolution due to greater numerical dispersion. Numerical error can be reduced by applying finer model resolutions at the cost of increased computation time. Table 2

Simulations with the Dynamic Wave Model
The previous example showed that the KWM gives unrealistic predictions of a flood wave routed through a prismatic channel with mild slope. This is largely due to the exclusion of inertial forces by truncating the momentum equation. Unsteady flow models that solve complete forms of Eq. (5), known as Dynamic Wave Models (DWMs), can provide more realistic simulations of unsteady flow.
The rivr package provides a DWM engine accessed in a similar manner to the KWM engine. The DWM engine uses the MacCormack predictor-corrector scheme (MacCormack, 2003) to provide solutions that are second-order accurate in space. As shown in Figure 7, flood waves modeled with the DWM attenuate as they travel through the channel.
The DWM engine differs from the KWM engine in that the downstream boundary condition must be known. The rivr package employs the Method of Characteristics (MOC) to resolve fluid flow across the upstream and downstream boundaries. By using MOC, the DWM engine can accommodate a variety of boundary specifications; the upstream and downstream boundaries can be individually specified in terms of either flows or water depths. The downstream boundary can alternatively be specified as a zero-gradient condition, which allows flow to be routed out of the channel by assuming that the flow or water depth is constant across the final (downstream) two nodes in the model domain. This results in "smearing" of the flood wave at the downstream boundary, but in some practical cases this may be preferable to direct specification of depth or flow. When depths are specified at the boundaries, the characteristic equations are solved directly; when flows are specified, a Newton-Raphson scheme is used to iteratively solve for depths and velocities simultaneously.
route_wave(slope, mannings, Cm, g, width, sideslope, initflow, boundary, downstream, timestep = dt, spacestep = dx, numnodes = numnodes, monitor.nodes = mpoints, monitor.times = mtimes, engine = "Dynamic", boundary.type = "QQ") The DWM engine, like the KWM engine, uses an explicit scheme and is subject to stability constraints. In fact, the DWM engine is more prone to numerical instability compared to the KWM engine and typically requires Courant numbers on the order of 0.06 (i.e. an order-of-magnitude increase in the required temporal resolution compared to the KWM engine for a given spatial resolution). This disadvantage is offset by a significant reduction in numerical error provided by the advanced finitedifferencing scheme. As shown in Figure 8, spatial resolutions of 1,000 feet provide accurate predictions at 50,000 feet downstream, at the cost of severe instability near the downstream boundary. Resolutions of 500 feet provide stable solutions throughout the full channel extent for a modest computation cost, as shown in Table 3.

Simulations with discontinuous boundary conditions
The DWM engine can accommodate boundary conditions that specify zero flows (velocities) but nonzero depths, such as can occur by the sudden closure of sluice gates. The last unsteady-flow example illustrates use of the DWM engine to simulate unsteady flow conditions following sudden closure of a sluice gate at the downstream end of a trapezoidal spillway. This example mimics example 14-1 qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q in Chaudhry (2007). The upstream boundary condition is specified as a constant depth to represent water level in a large reservoir. Example code for specifying discontinuous boundary conditions is shown below.

Opportunities for extension
The rivr package was designed such that additional analysis capabilities could be implemented in future versions with minimal modification to the existing code base. A wide variety of extensions could be supported by the existing interface, both in terms of numerical methods and problem complexity.
In particular, there are many avenues for contribution to the unsteady flow analysis capabilities of rivr.
Additional numerical schemes for unsteady flow analysis in prismatic channels-such as implicit finite-difference schemes and finite-element methods-could be implemented with only superficial changes to the current function interface provided by route_wave, i.e. new algorithms could be made accessible through additional keyword options in existing arguments provided the underlying engines followed the input template used by the KWM and DWM engines. Variable-step methods could also be implemented with minimal changes, most simply by allowing the user to define an initial or "default" time step and using an interpolation function to compute the upstream boundary condition at intermediate time steps when required. The authors consider the degree of effort required to implement such additional numerical schemes in rivr to be on par with typical expectations for an individual project assignment in a graduate-level computational hydraulics or numerical methods course.
Support for prismatic channels with arbitrary or compound geometry could be accomplished with only minor restructuring of the channel geometry function, but would introduce additional complexity to the specification of channel roughness (Yen, 2002) and the computation of normal and critical depth (Chaudhry and Bhallamudi, 1988;Younis et al., 2009). Support for compound channels is currently being considered for the next major version release of rivr.
Specification of variable channel geometry and slope, additional friction terms (e.g. bedform drag, hydraulic structures) and additional source or sink terms (e.g. hillslope runoff, infiltration) would require some modification to the current implementation but would largely consist of allowing users to specify channel properties as vectors rather than as single values, and modifying the underlying C++ code to access the appropriate vector elements when computing solutions at individual nodes. These modifications would best be developed prior to, or in conjunction with, the implementation of an implicit scheme. Such changes would, however, add considerable complexity to the formulation of unsteady flow problems and could potentially distract from the primary purpose of the package as an accessible tool for teaching. Support for channel networks would likely require the development of a new function interface that would essentially couple three or more instances of route_wave (e.g. two channels merging into one) and provide a keyword argument for selecting different methods for resolving flow depth and velocity at the junctions. Implementing such functionality would likely be a much larger endeavor compared to the other extensions discussed, and could pose additional challenges in terms of obtaining reasonably-fast solutions for moderately-large channel networks. It is important to note, however, that discussion of channel network modeling is often neglected in open-channel hydraulics courses due to the high complexity of problem formulation and solution algorithms; making channel network solution algorithms demonstrable through rivr would therefore be of considerable educational value.

Concluding remarks
This article describes rivr, a first-of-its-kind R/C++ package for one-dimensional open-channel flow computation with educational applications to hydraulic engineering degree programs and related disciplines. The package provides a reliable and flexible toolset for modeling open-channel flows via a simple function interface. Model outputs are formatted to facilitate analysis, tabulation, visualization and export. Example computations of backwater curves and flood wave routing are provided to demonstrate functionality and validate model results, and potential extensions to the package are discussed. This article, along with the package documentation, provides information on the governing equations and numerical implementations as well as a suite of use examples. The R Journal Vol. 7/2, December 2015 ISSN 2073-4859