Research article 18 Aug 2020
Research article  18 Aug 2020
A new form of the SaintVenant equations for variable topography
 ^{1}National Center for Infrastructure Modeling and Management, The University of Texas at Austin, Austin, Texas, USA
 ^{2}Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA
 ^{1}National Center for Infrastructure Modeling and Management, The University of Texas at Austin, Austin, Texas, USA
 ^{2}Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA
Correspondence: Ben R. Hodges (hodges@utexas.edu)
Hide author detailsCorrespondence: Ben R. Hodges (hodges@utexas.edu)
The solution stability of river models using the onedimensional (1D) SaintVenant equations can be easily undermined when source terms in the discrete equations do not satisfy the Lipschitz smoothness condition for partial differential equations. Although instability issues have been previously noted, they are typically treated as model implementation issues rather than as underlying problems associated with the form of the governing equations. This study proposes a new reference slope form of the SaintVenant equations to ensure smooth slope source terms and eliminate one source of potential numerical oscillations. It is shown that a simple algebraic transformation of channel geometry provides a smooth reference slope while preserving the correct crosssection flow area and the total Piezometric pressure gradient that drives the flow. The reference slope method ensures the slope source term in the governing equations is Lipschitz continuous while maintaining all the underlying complexity of the realworld geometry. The validity of the mathematical concept is demonstrated with the opensource Simulation Program for River Networks (SPRNT) model in a series of artificial test cases and a simulation of a small urban creek. Validation comparisons are made with analytical solutions and the Hydrologic Engineering Center's River Analysis System (HECRAS) model. The new method reduces numerical oscillations and instabilities without requiring ad hoc smoothing algorithms.
The SaintVenant equations (SVEs) for onedimensional (1D) river modeling are typically presented with pressure forcing terms of either (i) gradients of the water surface elevation or (ii) thalweg bottom slope combined with gradients of the water depth. In this study, we demonstrate a new form using a reference slope (S_{R}) and its associated depth (h_{a}), which are shown to be algebraically identical to the two standard forms of the SVEs. The new forms provide greater flexibility in addressing numerical convergence issues associated with modeling discontinuous bottom slopes. A key point of this paper is that precise representation of the thalweg bottom slope (S_{0}) and hydrostatic pressure gradients ($\partial {h}_{\mathrm{0}}/\partial x$) is not necessary to correctly represent variable topography. Indeed, the splitting point for representing the forcing Piezometric pressure gradient as a body force (defined by a slope) and a residual head gradient term is free choice in a simple algebraic substitution. Different choices for the splitting provide different body force directions and lead to different forms of the SVEs – all of which are valid representations of variable topography and do not constitute a smoothing of topography. We will show that it is possible to use a smooth slope (body force) term in the SVEs without actually smoothing the topography. Herein, this smooth slope term will be designated as the reference slope, S_{R}, to distinguish it from the traditional thalweg bottom slope, S_{0}.
The two common differential forms of the SVEs are
and
where Q is the flow rate, A is the crosssection area, η is the water surface elevation, h_{0} is the thalweg depth, S_{0} is the thalweg bottom slope, and S_{f} is the friction slope that represents the local energy gradient. Equation (1) can be envisioned as using the Piezometric head gradient to force the flow, as shown in Fig. 1a. In contrast, Eq. (2) can be envisioned as splitting the Piezometric gradient into a body force in the bottom slope direction and a hydrostatic head gradient, as shown in Fig. 1b. Note that both equations are valid for variable topography despite only the second equation explicitly representing the bottom slope and thalweg depth hydrostatic pressure gradients.
The two SVE forms of Eqs. (1) and (2) are algebraically identical using the identity
Heuristically, we can propose a more general identity of
where S_{R} is an arbitrary reference slope and h_{a} is an associated depth that will be defined in Sect. 3.2 below. For an introductory exposition, $\partial {h}_{\mathrm{a}}/\partial x$ is merely the residual implied for a given $\partial \mathit{\eta}/\partial x$ and arbitrary S_{R}. Applying Eq. (4) to (1) provides the following:
where the terms gAS_{R} and $gA\frac{\partial {h}_{\mathrm{a}}}{\partial x}$ are algebraically equivalent to those in Eqs. (1) and (2). Clearly, if we let S_{R} equal S_{0}, then h_{a} equals h_{0} and we recover Eq. (2). Furthermore, if we let S_{R} equal 0, then h_{a} equals η and we recover Eq. (1). The equations are identical with these substitutions, so it follows that using a reference slope of zero (S_{R}=0) must exactly represent the same topographic variability as using a reference slope that mimics the topographic slope (S_{R}=S_{0}) as long as h_{a} is correctly defined consistent with Eq. (4). That is to say that from simple algebra the use of the real S_{0} in the SVEs is not required to capture effects of topographic variability as long as the depth gradient term is correctly redefined as something other than the thalweg depth gradient and the crosssection flow areas are correctly computed.
From the arguments above, the effects of varying bottom topography are captured by S_{R}=0 and h_{a}=η, which implies we are also free to introduce any other (preferably smooth) S_{R} into Eq. (5) without altering the underlying representation of variable topography. An example is illustrated in Fig. 2. As the splitting defined in Eq. (4) makes Eq. (5) algebraically identical to Eqs. (1) and (2), the introduction of a smooth S_{R} does not reflect “smoothing” of the topography. It merely reflects a decision on whether effects of nonsmoothness will reside solely in solution variables A and h_{a} or will also be forced as a nonsmooth source term in S_{0}.
We would like to use an a priori smooth S_{R} in a computational model rather than the actual thalweg S_{0} because of what happens to S_{0}(x) and $\partial {h}_{\mathrm{0}}\left(x\right)/\partial x$ for topography varying sharply over short distances, as illustrated in Fig. 3. From a physics perspective, using S_{0} to split the Piezometric head is an intuitive way to describe the local interplay of pressure with the bottom slope. Furthermore, S_{0} has the advantage of being readily reduced to a kinematic wave equation where S_{f}=S_{0}, which has some advantage in multipurpose codes. However, from a numerical modeling perspective, using S_{0} has a significant limitation based on its smoothness. If the water surface is smooth, then nonsmooth S_{0}(x) requires the numerical solver to produce a compensating nonsmooth h_{0}(x), i.e., requiring a wellbalanced scheme (see Sect. 2). If we can discard our (wrong) intuition that the S_{0} form must somehow better represent sharply variable topography – i.e., recognizing the algebraic equivalence of Eq. (5) with Eqs. (1) and (2) – it follows that the splitting of the Piezometric head to include a body force that is everywhere exactly aligned with a sharply varying S_{0} is (from a numerical perspective) merely creating unnecessary complexity in the governing equation source term that requires compensating complexity in the solution algorithm. In contrast, by requiring S_{R} to be smooth, we can ensure the h_{a} solution is also smooth for a smoothly varying free surface.
The use of S_{R} rather than S_{0} in the governing equations can perhaps be better understood if we think of the slope in Eq. (4) as representing simply a portion of the overall Piezometric pressure gradient that can be extracted from $\partial \mathit{\eta}/\partial x$ and treated as a body force that varies gradually along the channel. Hence, in Fig. 3 we are not interested in separating out the details of the sharply varying slope changes in the local topography but instead prefer a body force term that aligns with the mean slope over some larger spatial scale, e.g., as in Fig. 2.
In this paper, we examine the effect of variable crosssection geometry on numerical solutions of the SVEs and propose a new reference slope approach that can be inferred from the above arguments. The use of a reference slope as a body force direction to split the Piezometric head gradient term ensures (i) the slope used in the discrete source term is smooth, (ii) the variable geometry is correctly retained, (iii) the fundamental governing equations are preserved, and (iv) an SVE numerical algorithm developed using S_{0} is essentially unchanged. We further demonstrate that bottom slope discontinuities are a cause of problems in finitedifference forms of 1D SaintVenant equations with subcritical flow. Of course, this idea will not be a surprise to many modelers who routinely remove troublesome cross sections or smooth their topography; however, the concept does not appear to have been conclusively demonstrated in the literature. More importantly, we show that the problem is inherent in the traditional formulation of the governing equations using the thalweg bottom slope, S_{0}, which is usually computed as the slope between the lowest points in two adjacent river cross sections. Problems associated with slope discontinuities can be fixed within the governing equations by the careful selection of smoothly changing reference elevations along the channel, z_{R}(x), which results in smooth reference slopes, S_{R}(x), and the redefinition of the thalweg depth (h_{0}) as a depth associated (h_{a}) with the reference elevation. Note that h_{a} is not necessarily any characteristic depth of the flow and is only indirectly related to the hydrostatic pressure.
The approach proposed herein can be implemented within any SaintVenant model as it is entirely independent of the solution algorithm; however, implementation does require rewriting code for the relationships between crosssection area, wetted perimeter, and the depth variable of the solution. Note that the new approach includes a redefinition of the depth variable (traditionally the maximum depth, h_{0}) as the depth associated (h_{a}) with the reference elevation, z_{R}.
Onedimensional (1D) hydrodynamic models using the SaintVenant equations (SVEs) are widely employed for studying both natural streams and manmade channels (e.g., MartinezAranda et al., 2019; Sanders, 2001). It is widely recognized that numerical solutions of the SVEs are prone to spurious oscillations in the freesurface elevation unless particular care is taken in the numerical formulation and/or the problem definition (e.g., Nujic, 1995; Tseng, 2004). Numerous techniques and special numerical schemes have been previously designed to overcome unwanted numerical oscillations caused by discontinuous geometries and boundary conditions (e.g., Zhou et al., 2001; Liang and Marche, 2009). These approaches typically rely on the concept of a wellbalanced discrete form (Greenberg and LeRoux, 1996) as discussed in a comprehensive review by Kesserwani (2013) and further elaborated by Hodges (2019).
Although wellbalanced schemes are relatively robust in handling discontinuous boundary conditions, they have not been extensively applied in water resource models to simulate the regionaltocontinentalscale river networks or stormwater systems for megacities. The rapidly varying and discontinuous S_{0} in natural systems can significantly increase the difficulty and computational burden of obtaining a wellbalanced method (Schippa and Pavan, 2008). Hence, when a largescale openchannel model develops oscillations and/or instabilities, practitioners may resort to the traditional approach of removing cross sections or smoothing bathymetry to mitigate oscillatory or unstable solution behavior (Tayfur et al., 1993). Such ad hoc efforts can be effective as they address a major cause of such oscillations and instabilities (discontinuous topography), but they inherently reduce the fidelity of the simulation.
Oscillations and instabilities can be induced in any numerical solution of a boundary initial value problem by the inclusion of nonsmooth source terms; i.e., if we consider an advection equation of the form
where σ is a nonhomogeneous source term, a fundamental theorem for differential equations provides that a unique solution cannot be guaranteed to exist unless the source term is Lipschitz continuous (e.g., Iserles, 1996). That numerical instabilities are often caused by nonsmooth source terms is not a new observation. A wide variety of numerical schemes have been developed to address this issue, including, e.g., extensive work on wetting/drying (Liang and Marche, 2009; Song et al., 2011), positivitypreserving methods for coupled models (Singh et al., 2015), and implicit schemes that address stiffness of the nonlinear friction term (Xia and Liang, 2018). The literature in this area is vast – particularly if both 1D and 2D models are considered. For the present purposes, we focus on only one part of the source term, S_{0}, whose nonsmoothness has previously been treated as a problem to be handled rather than as a problem that can be directly eliminated in the governing equations. Existing wellbalanced schemes (see reviews noted above) seek to compensate for nonsmoothness of all parts of the source term in the structure of the numerical discretization. Arguably, if the slope term is guaranteed smooth, then a wellbalanced scheme should be simpler to create. In general, when the thalweg bottom slope (S_{0}) appears as a source term in the SVEs, it should be a priori Lipschitz smooth or oscillations and instabilities should be expected. For natural systems, S_{0}(x) is typically defined using the maximum channel depth at each surveyed cross section, which is rarely a smooth function (unless the distance between cross sections is large compared to bottom elevation variability). When cross sections are surveyed at short distances, S_{0} will tend to have significant variability. It follows that the use of S_{0} has the undesirable property that smaller Δx (i.e., resolving a river with more detailed survey data) will increase the nonsmoothness in this source term of the momentum equation, resulting in a model that is unlikely to converge under a grid refinement test. It is not surprising that S_{0} smoothness, when it occurs in a model of a natural river channel, is typically the result of relatively long separations (Δx) between crosssection surveys that ensures that the discrete d^{2}z_{0}∕dx^{2} is small. Thus, removing cross sections can be an effective mitigation technique because it increases Δx and effectively smooths S_{0}. In general, models discretized with higherresolution river surveys (smaller Δx) will have greater nonsmoothness in S_{0} and develop more oscillation and instability issues. In essence, our models get worse as our boundary condition data get better.
The problems associated with S_{0} can be understood by considering the identity in Eq. (3) for a channel with subcritical flow in which the freesurface curvature is expected to be negligible, i.e., ${\partial}^{\mathrm{2}}\mathit{\eta}/\partial {x}^{\mathrm{2}}\sim \mathit{\u03f5}$, where $\mathit{\u03f5}\ll \text{d}{S}_{\mathrm{0}}/\text{d}x,{\partial}^{\mathrm{2}}\mathit{\eta}/{\partial}^{\mathrm{2}}x$. Taking the alongchannel gradient of Eq. (3) implies that
Thus, forcing with a nonsmooth dS_{0}(x)∕dx will require nonnegligible curvature of the response variable, h_{0}(x), whose gradient is also a forcing function of the nonlinear equation. Feedback can easily build and cause successive overshoot/undershoot effects, producing oscillations and nonconvergence in a nonlinear solver. In contrast, Eq. (4) can be invoked with dS_{R}∕dx guaranteed to be small, which implies that ${\partial}^{\mathrm{2}}{h}_{\mathrm{a}}/\partial {x}^{\mathrm{2}}$ will also be small even when $\partial A/\partial x$ is nonsmooth. As a practical matter, any S_{0} with a discontinuous discrete first derivative (i.e., discontinuities in the second derivative of the thalweg elevation, d^{2}z_{0}∕dx^{2}) will be Lipschitz discontinuous and should not be directly discretized in an SVE solution with Eq. (2). Although approximate numerical solutions of equations with nonsmooth S_{0}(x) can sometimes be attained for models with sufficient damping, such solutions are questionable as they do not have rigorous mathematical foundations.
Arguably, nonsmoothness in S_{0} can be handled in one of four ways: (i) smoothing the geometry – hence solving for flows that do not match the real system; (ii) applying ad hoc smoothing within the flow solution – i.e., adjusting the physics to remove numerical instabilities; (iii) adjusting the numerical discretization scheme to compensate for nonsmoothness – e.g., the wellbalanced concept; or (iv) adjusting the governing equations to ensure that any slope in the source term is smooth without modifying the solution physics, the channel geometry, or the numerical discretization scheme. It should be obvious that Eq. (1) is the extreme example of the last approach – replacing S_{0} and $\partial {h}_{\mathrm{0}}/\partial x$ with $\partial \mathit{\eta}/\partial x$ ensures that S_{0} does not occur in the governing equations and cannot destabilize the solution. To our knowledge, the last approach (as used herein and illustrated in Fig. 2) has not been previously proposed or analyzed in the literature. Nevertheless, as shown below, it provides a simple method that can be readily adapted into existing hydrodynamic models.
For brevity, we will limit our focus herein to subcritical flows – backwater tends to smooth the effects of slope discontinuities and thus we expect smooth solutions for flow rate and freesurface elevation despite nonsmooth geometry. Nevertheless, common SVE solvers can exhibit oscillatory, nonconvergent behavior even in simple subcritical flows when geometry is not smooth. In the following, it will be obvious that the mathematical theory applies directly to supercritical and transcritical flows as well, but evaluating model performance under the breadth of possible transcritical conditions (including nonsmooth jumps) necessarily requires more analyses than is practical in a single paper.
3.1 SPRNT
The Simulation Program for River Networks (SPRNT) code for unsteady SVE river networks is used and modified herein. The baseline for this code models momentum using Eq. (2), which is coupled to the solution of continuity,
where q_{ℓ} is a lateral inflow per unit length. Note that significantly nonsmooth q_{ℓ}(x,t) can provide another source of numerical oscillations and instability (Kuiry et al., 2010). As the main focus of this study is the slope source term in the momentum equation, the lateral inflows and their effects are neglected by setting q_{ℓ}=0 everywhere.
In the SPRNT momentum equation, the friction slope is represented using the Chezy–Manning form as
where P_{w} is the wetted perimeter of a cross section and n is the Gauckler–Manning–Strickler roughness. Although, there are other methods for treating frictional losses (e.g., Decoene et al., 2009; Burguete et al., 2007), the Chezy–Manning form remains popular due to its simplicity.
The baseline model uses the thalweg elevation (z_{0}), the thalweg depth (h_{0}), and the thalweg bottom slope (S_{0}) as
and
SPRNT is an opensource, 1D hydrodynamic solver using the fully implicit Preissmann numerical scheme (Preissmann, 1961) with the Newton–Raphson iteration and computational acceleration techniques developed from very largescale integration (VLSI) semiconductor design. Details on the baseline SPRNT model and its application to large river networks are provided in Liu and Hodges (2014) and Yu et al. (2017).
3.2 Reference slope method
We introduce a new reference slope (RS) method through a transformation and redefinition of geometry in the SaintVenant equations, as discussed in Sect. 1. In place of the conventional h_{0} and z_{0}, we define a reference elevation (z_{R}) and its associated depth (h_{a}) as shown in Fig. 4. These provide a relationship with the freesurface elevation (η) defined as
Note that z_{R} is arbitrary, so h_{a} may be either greater than or less than the thalweg depth, h_{0}, at a given location. As shown in Fig. 4, it is convenient to define the reference height, h_{R}, in relation to z_{R} and the true bottom elevation, z_{0}, by
Thus, the conventional h_{0} and z_{0} are recovered with
and
We define the reference slope (S_{R}) as the downstream slope of z_{R}:
Using Eqs. (12) and (16) in Eq. (2) provides Eq. (5), which is more conveniently written as
The above is identical to Eq. (2) with the simple substitution of h_{a} and S_{R} for h_{0} and S_{0}. In this formulation, the definition of z_{R}(x) is arbitrary, so we can a priori require a definition such that S_{R}(x) is smooth. A trivial choice that is guaranteed smooth is z_{R}(x)=constant, which returns S_{R}(x)=0 and the SaintVenant equations in the form of Eq. (1). However, this form with $\partial \mathit{\eta}/\partial x$ for the entire pressure term is known to cause numerical stiffness issues for large ranges in η, e.g., the elevation change in a river from its mountain source to a coastal plain (Liu and Hodges, 2014). Using the conventional S_{0} in Eq. (2) reduces this problem as the range of h_{0} is inherently confined to the local water depths rather than the underlying topography. In the RS method, the range of h_{a} is tied to the range of water depths and the selection of z_{R}; thus, for present purposes, we are interested in nontrivial definitions of z_{R} that (i) are close to z_{0} to maintain a small range of h_{a} values and (ii) provide smooth S_{R}(x). An arbitrary z_{R} that is far from z_{0} or nonsmooth is of little interest as it holds no theoretical or practical advantage over the Eq. (1) approach implied by S_{R}(x)=0.
If S_{R}(x) is required to be smooth, then the source term of the equation can be guaranteed smooth as long as S_{f}(x) is smooth, which is typically true as long as the solution, Q(x), is smooth. Note that in extreme cases of geometric discontinuity the combined values of n, P_{w}, and A in Eq. (9) can cause a nonLipschitz friction term; thus, the RS method cannot guarantee that the entire source term is smooth. Numerical solution methods are usually robust against discontinuities in n and P_{w} as they are coefficients of the solution variables {A,Q}. More subtle problems might arise due to discontinuities developed in the ${Q}^{\mathrm{2}}{A}^{\mathrm{10}/\mathrm{3}}$ ratio in Eq. (9); countering incipient instabilities from this term requires other numerical strategies (e.g., Xia and Liang, 2018).
A critical change required by the introduction of h_{a} is that the conventional geometric auxiliary relationships of A=f(h_{0}) and P_{w}=f(h_{0}) must be transformed into A=f(h_{a}) and P_{w}=f(h_{a}). That is, once we change the depth in our equation from h_{0} to h_{a}, we must reindex the geometry. In general, for known functions A(h_{0}) and P_{w}(h_{0}) this is a trivial transformation:
and
However, the implementation in an existing code is not necessarily as simple as the above equations suggest. For example, in Fig. 4, the condition A=0 occurs when h_{a}=h_{R}, i.e., a nonzero value as compared to h_{0}=0 with conventional geometry. Unfortunately, model developers typically have ad hoc wetting/drying treatments that are introduced as h_{0}→0 or for h_{0}<0. Such treatments need to be modified to be deployed as h_{a}→h_{R}, which introduces the added complication that h_{R} is negative when z_{R}>z_{0}. Note that the new geometry does not require altering the wetting/drying algorithm itself or, for that matter, any other solution algorithm – only the actual geometry definitions require alteration. These relatively straightforward changes can be contrasted with the effort needed to provide a wellbalanced numerical discretization scheme for the conventional S_{0} representation of geometry (e.g., Kesserwani, 2013).
The modification of the SPRNT code to implement the above RS method will be known as SPRNTRS. The SPRNT and SPRNTRS source codes are available in an opensource repository (Liu, 2014). Note that the solution algorithm for SPRNTRS is identical to that of SPRNT; the only code changes are for the new geometry definitions for h_{a} and S_{R} that replace h_{0} and S_{0} in the original algorithm. This simple geometry replacement strategy is effective because Eq. (17) is identical to Eq. (2) except for the change in nomenclature to h_{a} and S_{R}.
3.3 Generating a smooth S_{R}(x)
The z_{R}(x) and hence S_{R}(x) are arbitrary choices in the RS method but should be generated for the smoothness of S_{R}(x) along the reach. In synthetic reach test cases and analytical test cases (described below), the channels are a priori either specified with a uniform S_{R} or produced by a different order of splines. For our urban creek test case, the z_{R} is generated with an approximating cubic Bspline (de Boor, 2001) based on thalweg, z_{0}(x), elevations. There are a variety of possible ways to generate smooth S_{R}(x), but applying approximating cubic Bsplines to the z_{0}(x) is convenient because the slope is guaranteed to be locally smooth as long as the knot spacing in the Bspline is everywhere larger than the spacing between cross sections. It should be emphasized that an exact spline fitting of all the thalweg data (i.e., knots at all the cross sections) will be smooth at scales finer than the crosssection spacing but nonsmooth at the model's discretization scale. That is to say that the exact cubic spline fitting of z_{0}(x) does not reduce discontinuities at the discretization scale – only an approximate fitting associated with coarser scales than the crosssection spacing will be effective. It follows that there is some (limited) choice in the selection of the subset of z_{0}(x) used as the spline knots, with different sets producing slightly different {z_{R}(x),h_{R}(x)} over the domain. Each set is algebraically identical to the underling geometry, so the generated solutions should be identical within the machine truncation error as long as the z_{R}(x) are sufficiently smooth. Implications of the method chosen for generating z_{R}(x) are discussed in Sect. 5.3. Further details and test cases are provided in Yu et al. (2019b).
3.4 HECRAS for model validation
The baseline SPRNT has been previously shown to have excellent agreement with the Hydrologic Engineering Center's River Analysis System (version 5.0.7) – known as HECRAS. Liu and Hodges (2014) showed SPRNT simulations agreed with HECRAS with ≤3 % difference in water depth solution when using both prismatic cross sections and nonuniform channels. Thus, HECRAS provides a reasonable model for testing and validating SPRINTRS. We would have preferred to use a single model with and without the RS method for such model–model comparisons; however, HECRAS is a closedsource proprietary model, so we could not directly implement and test the RS method in that code. Conversely, as expected by the discussions in Sect. 2, the baseline SPRNT model is oscillatory and nonconvergent on the highly discontinuous geometry of our test cases due to its use of the S_{0} approach, so it cannot be directly used for before and after comparisons of RS. Thus, simulations using SPRNTRS are compared to HECRAS simulations for validation and insight.
HECRAS provides a convenient validation model for three reasons. Firstly, it is a widely accepted engineering model for riverreach simulations, (e.g., Wang et al., 2012; Giustarini et al., 2011; Aggett and Wilson, 2009). Secondly, it has been used as a validation model in numerous prior studies (e.g., Gichamo et al., 2012; Mejia and Reed, 2011; Horritt and Bates, 2002). Finally, unsteady HECRAS employs $\partial \mathit{\eta}/\partial x$ as the piezometric gradient rather than using $\partial {h}_{\mathrm{0}}/\partial x$ and S_{0}, which is one of the reasons it is relatively robust for nonsmooth geometry such as that used herein.
The performance of the RS method is demonstrated below through (i) a comparison to six analytical test cases from MacDonald et al. (1995) with Lipschitz continuous geometry, various prismatic shapes, and different formulations of S_{R}, (ii) seven synthetic test cases using Lipschitz discontinuous geometry, and (iii) an urban creek with complex crosssection geometry derived from physical surveys that include discontinuities anorderofmagnitude greater than those in the synthetic test cases.
3.5 Test cases – analytical solutions
Analytical solutions of six test cases with different channel shapes and bed slope formations from MacDonald et al. (1995) are used to show that SPRNTRS reproduces the correct water surface elevation regardless of the selection of S_{R}. These test cases are representative of the more comprehensive analysis provided in Yu et al. (2019a). The configuration details for each case are provided in Table 1, for which we adopt the nomenclature of MacDonald et al. (1995) for ease of comparison. The selected test cases have Lipschitz smooth geometric features that are represented in RS tests using both uniform and splined reference beds, as shown in Fig. 5. To illustrate the adaptability of the RS method, the UR3 and UT2 cases use splines that produce z_{R} very close to (but not identical to) the actual bed, whereas the other cases use uniform S_{R} or splines with greater differences. The uniform S_{R} in cases UR1, UT1, and VR1 is set to the average slope in each domain. With reference to Fig. 4, the differences between the channel bottom (z_{0}) and the reference bottom (z_{R}) shown in Fig. 5 imply channel bottom offsets (h_{R}) of varying complexity for the RS method, as shown in Fig. 6. The VR1 and VR2 cases also provide smooth changes in the channel width, which are shown in Fig. 7. The node spacing for all these tests is a uniform Δx=10 m. The boundary conditions follow MacDonald et al. (1995).
3.6 Test cases – synthetic channel reach
The analytical test cases above are designed to show that the RS method does not introduce approximations that affect the smooth solution. However, the true power of the RS method is in solutions of nonsmooth bathymetry in which most models using h_{0} and S_{0} have difficulty converging. To illustrate this aspect, we use a simple river reach with a randomly perturbed (discontinuous) bathymetry at various scales. As there are no analytical solutions for these tests, we use the HECRAS model for comparison. The simulations use timeinvariant boundary conditions with geometry defined by trapezoidal cross sections of uniform side slope, as detailed in Table 2. The channel bed offset (h_{R}) and thalweg slope (S_{0}) are illustrated in Fig. 8. The flow boundary conditions provide for a mild slope with a water surface profile that can be classified as an M1 gradually varying flow.
Case 1 is the baseline smooth channel with a uniform slope over the entire reach length. Cases 2 through 5 have a synthetic geometry developed by random perturbations of the bottom elevation of the baseline reach. Cases A and B have the identical smooth geometry to Case 1 but use different reference slopes for RS tests. The synthetic channel test reach is 1.58 km in length discretized into 80 uniform computational nodes with 79 channel segments (20 m per segment). The trapezoidal cross sections each have a 10.0 m bottom width and 63.4^{∘} sidewall slopes. Bottom roughness is fixed by a Manning's n value of 0.04 for all segments.
To develop the random perturbations of the bottom in the synthetic test cases, we begin with Case 1 (baseline) using a uniform ${S}_{\mathrm{0}}^{\left[\mathrm{1}\right]}=\mathrm{0.008}$ over the entire reach length. Here the superscript in square brackets indicates the case identifier. The set of bottom elevations for Case 1 are ${z}_{\mathrm{0}}^{\left[\mathrm{1}\right]}\left(x\right)$, which are smooth and linearly decreasing over the reach length. Cases 2–5 are similar channels with perturbed bottom elevations set by
where H(x) is a set of randomgenerated numbers within the range of $\mathrm{0.126}\le H\left(x\right)\le \mathrm{0.183}$. The upper and lower limits of H(x) were selected to prevent the occurrence of a locally adverse slope – such conditions can be handled by SPRNTRS but can cause convergence problems for some models. The α^{[c]} is a magnitude to generate a range of bottom displacements with ${\mathit{\alpha}}^{\left[c\right]}\in \left\{\mathrm{0.01},\mathrm{0.1},\mathrm{0.5},\mathrm{1.0}\right\}$ for $c\in \left\{\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}\right\}$, respectively. Cases 2–5 set the reference bottom elevations exactly equal to the baseline Case 1 physical bottom elevations; i.e.,
Thus, the SPRNTRS simulations for cases 2–5 use uniform S_{R} over the reach such that the bed offset (h_{R}) represents the physical geometric perturbations. Noting from Eq. (13) that h_{R} is the difference between the physical bottom (${z}_{\mathrm{0}}^{\left[c\right]}$) and the reference bottom (${z}_{\mathrm{R}}^{\left[c\right]}$), substituting the above relationships gives
For synthetic test cases A and B in Table 2, the actual channel bottom slopes are set to uniform values equivalent to Case 1; that is ${S}_{\mathrm{0}}^{\left[A\right]}={S}_{\mathrm{0}}^{\left[B\right]}={S}_{\mathrm{0}}^{\left[\mathrm{1}\right]}=\mathrm{0.008}$ with identical thalweg elevations of ${z}_{\mathrm{0}}^{\left[\mathrm{1}\right]}\left(x\right)$. However, the reference slopes for these cases are set to smaller and greater uniform values: ${S}_{\mathrm{R}}^{\left[A\right]}=\mathrm{0.004}$ and ${S}_{\mathrm{R}}^{\left[B\right]}=\mathrm{0.010}$. These two cases demonstrate that the RS method generates the same numerical solution as baseline Case 1 (solved at S_{0}) when S_{R} is set to an arbitrary value.
Forcing for all seven test cases is a constant inflow boundary of 283 m^{3} s^{−1} applied at the furthest upstream node. The downstream boundary condition is 5.0 m depth, which is a subcritical flow based on a normal depth of 4.95 m for S_{0}=0.008 and the inflow rate. Because a subcritical boundary condition allows upstream wave reflections, the downstream boundary was enforced at the end of a 180 m (9 node) buffer domain, which was adequate for reducing upstream wave propagation in unsteady flow solutions. Simulation results are reported after the models have reached a steady state and all oscillations associated with the initial conditions have dissipated. Solutions for the buffer segment are not included in the analyses below.
3.7 Test case – Waller Creek study site
The main stem of Waller Creek in Austin, Texas, USA, is used to examine the performance of the RS method for more complex conditions. The main stem of the creek drains an urban watershed of 14.3 km^{2} with total length of 10.7 km for the area illustrated in Fig. 9. Bathymetric survey data are available courtesy of the city of Austin (Fig. 10). The bathymetric data set includes 327 surveyed cross sections with spacing intervals ranging from 2.5 up to 178 m (mean of 33.5 m). The Manning's n of the channel (based on the city of Austin computations) varies from 0.02 to 0.06 throughout the system. The SPRNT model, in both its original and RS form, has numerical stability issues associated with close crosssection spacing. Arguably, these issues are related to sharp changes in A that lead to nonsmooth source terms despite the RS method; however, this issue requires further investigation. Similar numerical instability behavior can also be found in the HECRAS unsteady model and also causes divergent solutions. For the present work, we discarded 36 cross sections (11 % of the data set) that were closer than 10 m and merged these short reach lengths with the adjacent sections. An additional three cross sections were discarded and some channel roughness values were modified as they caused numerical instability in the HECRAS unsteady simulation (see Appendix for details). The resulting data set is 288 cross sections with spacing ranging from 10.1 to 184.9 m. The mean crosssection spacing is 37.2 m with a total reach length of 10.7 km. To limit our focus to subcritical flow, our analyses consider only the upper 8.3 km of the main reach (210 out of 288 cross sections), which eliminates a series of steppool transcritical elements in the downstream channel in which the HECRAS solution is strongly influenced by the ad hoc local partial inertia (LPI) algorithm (Fread et al., 1996; Brunner, 2016b). The smoothing introduced by LPI makes it difficult to draw conclusions from a comparison between SPRNTRS and HECRAS across transcritical locations.
A timeinvariant upstream inflow boundary condition is set to 25 m^{3} s^{−1} at the headwater cross section. To minimize the influence of subcritical reflections from the upstream inflow boundary, the first 10 computational nodes at the upstream are not included in the results analysis. Lateral inflows are set to zero for all test cases. A 300 m buffer section is added downstream of the test domain to reduce the influence of reflections from the downstream boundary condition. This buffer section uses the same cross section as the final downstream section of the data set with a bed slope (S_{0}) of 0.0033 and Manning's n of 0.04. The buffer section has a normal depth of 0.76 m at the 25 m^{3} s^{−1} inflow rate. The downstream boundary condition at the end of the buffer section is 0.7 m depth, which is subcritical and implies an M2 gradually varying drawdown in the vicinity of the outflow. These geometry and boundary conditions are identically applied to both SPRNTRS and HECRAS models.
The thalweg elevation, z_{0}(x), and the reference elevation, z_{R}(x), of the RS method (as determined by the approximate spline fit described above) for Waller Creek are shown in Fig. 11a. The z_{0} and z_{R} are visually similar with the former being somewhat more noisy. The elevation data sets provide similar overall reach slopes (uppermost cross section to lowermost cross section) of 0.0074 and 0.0077, respectively. Note that the approximate Bspline technique for generating z_{R}(x) does not force the overall reach slope to be identical. Because z_{R}(x) is mathematically arbitrary, there is no need to force an exact match. Although z_{0}(x) and z_{R}(x) are similar in Fig. 11a, S_{0}(x) from the raw data is discontinuous and varies over a wide range (up to 4 times the reach slope), as illustrated in Fig. 11b. Note that S_{0}(x) also includes negative slopes (i.e., adverse gradient sections) which can cause convergence problems for some numerical solvers. In contrast, as shown in Fig. 11c, the approximate cubic Bspline used to generate S_{R}(x) from z_{R}(x) provides a reference slope that is everywhere smooth and positive and remains close to the overall reach slope of 0.008. The slope range and model nomenclature for the Waller Creek test cases are provided in Table 3. The Lipschitz smoothness of S_{0} versus S_{R} can be better understood by evaluating the gradient of the slope, i.e., the second derivative of z_{0} and z_{R}, as shown in Fig. 11d. The S_{0} formulation clearly lacks smoothness in the higher derivative.
3.8 Analysis methods
To evaluate the performance of the RS method relative to conventional formulations, four depthbased indicators are employed, as described below. For these definitions, the control (superscript [C]) is the MacDonald et al. (1995) solution for the analytical test case and HECRAS results for the synthetic channel and Waller Creek test cases. Note that the synthetic tests use unsteady HECRAS, whereas the Waller Creek study uses comparisons to both steady and unsteady versions of the model. The test case (superscript [T]) is always the SPRNTRS simulation. These measures can be considered error metrics for the comparison to the analytical solutions and difference metrics for model–model comparisons.

Normalized difference (ρ). A nondimensional index to describe the local difference in depth can be defined as
$$\begin{array}{}\text{(23)}& {\mathit{\rho}}^{[T:C]}\left(x\right)={\displaystyle \frac{{h}_{\mathrm{0}}^{\left[C\right]}\left(x\right){h}_{\mathrm{0}}^{\left[T\right]}\left(x\right)}{{h}_{\mathrm{0}}^{\left[C\right]}\left(x\right)}},\end{array}$$where ${h}_{\mathrm{0}}^{\left[C\right]}\left(x\right)$ and ${h}_{\mathrm{0}}^{\left[T\right]}\left(x\right)$ are the local depth solutions from the control and test case results after steadystate conditions are achieved. The normalization scale is the local depth of the control case. Note the denominator is nonzero in the synthetic test case setup because the flow setup is an M1 gradually varying flow.

Absolute mean normalized difference (ζ). The mean of the absolute value of ρ(x) over the domain provides an overall nondimensional indicator of the depth error:
$$\begin{array}{}\text{(24)}& {\mathit{\zeta}}^{[T:C]}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{x=\mathrm{1}}^{N}\mid {\mathit{\rho}}^{[T:C]}\left(x\right)\mid ,\end{array}$$where N is the total number of cross sections. We use the absolute value so that positive errors do not cancel negative errors, and ζ is a representative scale of the discrepancy between models.

Mean absolute error (MAE). The overall dimensional error is characterized as
$$\begin{array}{}\text{(25)}& \text{MAE}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{x=\mathrm{1}}^{N}\mid {h}_{\mathrm{0}}^{\left[C\right]}\left(x\right){h}_{\mathrm{0}}^{\left[T\right]}\left(x\right)\mid ,\end{array}$$and the nondimensional form of overall error is
$$\begin{array}{}\text{(26)}& \begin{array}{rl}\text{MAE}& \phantom{\rule{0.25em}{0ex}}\text{(nondimensional)}\\ & ={\displaystyle \frac{\mathrm{1}}{N}}\sum _{x=\mathrm{1}}^{N}\mid {\displaystyle \frac{{h}_{\mathrm{0}}^{\left[C\right]}\left(x\right){h}_{\mathrm{0}}^{\left[T\right]}\left(x\right)}{{h}_{\mathrm{0}}^{\left[C\right]}\left(x\right)}}\mid .\end{array}\end{array}$$ 
Root mean square error (RMSE). A standard dimensional measure of the squared error is
$$\begin{array}{}\text{(27)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{x=\mathrm{1}}^{n}\left({h}_{\mathrm{0}}^{\left[C\right]}\right(x){h}_{\mathrm{0}}^{\left[T\right]}(x){)}^{\mathrm{2}}}.\end{array}$$The nondimensional form of RMSE is computed by the following equation:
$$\begin{array}{}\text{(28)}& \phantom{\rule{0.33em}{0ex}}\begin{array}{rl}\text{RMSE}& \phantom{\rule{0.25em}{0ex}}\text{(nondimensional)}\\ & =\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{x=\mathrm{1}}^{n}({\displaystyle \frac{{h}_{\mathrm{0}}^{\left[C\right]}\left(x\right){h}_{\mathrm{0}}^{\left[T\right]}\left(x\right)}{{h}_{\mathrm{0}}^{\left[C\right]}\left(x\right)}}{)}^{\mathrm{2}}}.\end{array}\end{array}$$Both the MAE and RMSE are also reported in nondimensional form where the normalization scale is presented.
4.1 Analytical test cases
The water surface elevations for the analytical solutions and SPRNTRS simulations are shown in Fig. 12. Visually, the analytical and simulated results across all six cases are identical. Error metrics following Sect. 3.8 are provided in Table 1. The normalized differences (ρ) are less than 1 % and are consistent with absolute errors of O(10^{−3}) m, which are negligible compared with the water depths ≥1 m. The spatial distributions of the normalized error are shown in Fig. 13. By comparing this figure with Fig. 5, it can be seen that ρ(x) fluctuates with the change in bed slope. Similar behavior can also be found for model results reported in MacDonald et al. (1995).
4.2 Synthetic test cases
Results for the baseline synthetic test, Case 1, are shown in Fig. 14. The SPRNTRS method produces visually the same solution to HECRAS with z_{R}=z_{0}. Similarly, the comparison of model results for depth (h_{0}) for test cases 2–5 are visually indistinguishable, as shown in the left column of Fig. 15. The quantitative difference measures for the synthetic tests are provided in Table 5, and the spatial distributions of ρ(x) are illustrated in the right column of Fig. 15. Values for ρ(x) in cases 2 and 3 are slightly below zero (≈0.02 %) over the entire domain, indicating the SPRNTRS solution has a slightly higher water surface than the HECRAS solution for small perturbations in the bed slope. With the increased bottom perturbations in cases 4 and 5, the ρ(x) range is larger (and includes a positive range), but the bounding values are still trivial. The ζ and RMSE measures show that the nondimensional and dimensional overall differences are small. The MAE and RMSE climb slightly with the increasing h_{R} for cases 2 through 5 but remain below 3 mm. These depth RMSE values are negligible compared with the normal depth (4.95 m) of the baseline and well within reasonable truncation error differences for solvers using different numerical techniques. The model–model comparisons for test cases A and B also have trivial errors (Table 5), and further results are not shown as they are visually identical to those for baseline Case 1 illustrated in Fig. 14. Note that the RMSE for both cases is identical to Case 1, which indicates the solution for the SPRNTRS method with S_{R}≠S_{0} is very close to the baseline solution with S_{R}=S_{0}.
4.3 Waller Creek test case
Waller Creek has been simulated with SPRNTRS (denoted as WC_{RS} in the following figures), the HECRAS unsteady solver (WC_{HECU}), and the HECRAS steady solver (WC_{HECS}). Figure 16 shows water surface elevations for SPRNTRS and unsteady HECRAS. For clarity, the upper 40% of the domain (which has similar good behavior) is not shown. Figure 17 shows the spatial distribution of the normalized difference, ρ, for these simulations. The differences are roughly within ±4 % across the entire domain. The maximum and minimum difference both occur at two adjacent nodes close to 7800 m with 4.14 % and −3.07 %, respectively. Figure 18 provides a similar comparison of water surface elevations between SPRNTRS and the steady HECRAS case. The results are visually quite similar to the comparison with unsteady HECRAS. A direct comparison of surface elevations for unsteady and steady HECRAS does not provide any further insight and is omitted for brevity. However, to quantitatively evaluate the differences between SPRNTRS and HECRAS, it is useful to compute difference measures between the unsteady and steady HECRAS models themselves, as well as the differences between SPRNTRS and both models, as provided in Table 6. Overall, the SPRNTRS result have marginally better consistency with unsteady HECRAS than with steady HECRAS. Of greater importance is that the behavior of SPRNTRS relative to unsteady HECRAS has the same order of differences as the comparison of unsteady HECRAS to steady HECRAS. These results imply that the differences between SPRNTRS and unsteady HECRAS are reasonable for the different numerical methods given the geometric variability of Waller Creek.
5.1 Validation of the RS method
The RS method is a simple algebraic transformation of the governing equations, and the answer to the principal question “does it work” is implied by our inability use the baseline SPRNT model (with S_{0}) as a control model (see Sect. 3.4). Invariably, discontinuous topography for SPRNT without RS caused either an oscillatory solution or numerical instability. In contrast, both HECRAS (using η) and SPRNTRS (using S_{R}) provide stable, nonoscillatory solutions.
The analytical results in Sect. 4.1, supplemented by additional results in Yu et al. (2019a), validate the SPRNTRS method for the simulation of smoothly varying channel morphologies that are Lipschitz continuous at the discretization scale. We have experimented with both uniform and splined S_{R} for these tests. For both types of simulations, we observe errors relative to physical experiments that are comparable or smaller than those shown in the numerical validation studies of MacDonald et al. (1995). These results imply that the transformation from the conventional h_{0} and S_{0} form of the SVEs to the h_{a} and S_{R} form of Eq. (17) is a valid algebraic step that can be implemented in a numerical solver and is an alternative for representing smooth geometries.
The synthetic test cases in Sect. 4.2 serve two purposes. Firstly, cases A and B compared to baseline Case 1 show that the numerical solution does not depend on a particular choice of S_{R}. The arbitrary selection of an S_{R}≠S_{0} results in identical solutions to S_{R}=S_{0}. Secondly, the synthetic test cases show that the SPRNTRS method can be applied with nonsmooth geometry at the discretization scale (i.e., random perturbations of the physical bottom slope), which caused nonconvergent behavior in the baseline SPRNT model. As a control, we have compared SPRNTRS with the accepted HECRAS model that remains stable for these test cases as it is solved with the piezometric pressure gradient rather than splitting into S_{0} and the gradient of h_{0}. The results indicate that SPRNTRS provides numerical solutions that are nearly identical to HECRAS for the nonsmooth geometry test cases. Thus, using a Lipschitz smooth S_{R} provides a stable numerical solution for nonsmooth geometry without altering the physical representation of nonsmooth geometry.
The Waller Creek test case in Sect. 4.3 provides a more challenging comparison of SPRNTRS to HECRAS. For this test case, the geometry discontinuities include adverse slopes and local S_{0} that are ±400 % of the reachaverage slope, which contrasts with perturbations of ±30 % used in the synthetic test cases of Sect. 4.2. Again, SPRNTRS is shown to be close to the unsteady HECRAS solution. The model differences are within reasonable ranges, as illustrated by the fact that they are similar to the differences between HECRAS steady and unsteady versions. Nevertheless, it remains possible that the minor differences between HECRAS and SPRNTRS are caused by a latent defect in coding the RS method or in SPRNT itself, but it is difficult to envision how such a defect could occur without also appearing in the analytical and synthetic test cases. A simpler and more compelling explanation is with the linear approximations used in unsteady HECRAS that are not present in SPRNTRS. Specifically, Brunner (2016a) notes that for computational efficiency and to reduce “troublesome convergence problems at discontinuities in the river geometry,” the unsteady HECRAS code uses a linearization technique developed by Liggett (1975) and Chen (1973) – note that the latter document is cited by Brunner (2016a) but was not available to us. It seems likely that strong geometry discontinuities in the Waller Creek test case would be affected by this linearization, which arguably would lead to artificial smoothing of the water surface profile by HECRAS. Unfortunately, we do not have direct access to the HECRAS code and thus rely on the discussion of HECRAS stability in the literature (Hicks and Peacock, 2005; Sharkey, 2014) and the methodology in HECRAS manuals (Brunner, 2016a, b).
5.2 Why not just use $\partial \mathit{\eta}/\partial x$?
One might wonder whether S_{R} or S_{0} is at all necessary when we could clearly just retain $\partial \mathit{\eta}/\partial x$ in the SVEs rather than using any split form. To understand the value of S_{R}, it is worth considering why S_{0} is presently used. We have not been able to determine exactly when S_{0} was first used with the SVEs, but from a hydrology viewpoint S_{0} provides consistency between kinematic wave solutions (which use S_{0}=S_{f}) and the SVEs. Thus including S_{0} is a logical step when considering reducedphysics approaches (Di Baldassarre, 2012). Arguably, a wellchosen S_{R} that matches the largescale S_{0} will serve the same purpose. The S_{0} approach is also favored in models that are built on a “conservative” SVE form in which the hydrostatic pressure portion of the piezometric head gradient is abstracted into the advective gradient term (e.g., Sanders, 2001; Kesserwani, 2013). For these model, the advantage of the S_{0} form is that, when S_{0} equals 0 and S_{f} equals 0, the momentum equation can be written as a classic 1D homogenous advection equation, which is mathematically appealing. Our work in progress indicates that the S_{R} approach could be similarly adapted for a conservative form of the SVEs, but this issue remains speculative.
Although the utility and simplicity of the η approach is obvious, it has a key disadvantage when applied in largescale simulations. Over large distances the free surface is monotonically increasing upstream, which has consequences for employing implicit or semiimplicit numerical solutions in a continental river dynamics framework (Hodges, 2013). Briefly, when modeling a river system from an estuary (η∼0 m) to mountain headwaters (η∼10^{3} m), the solution variable, η, nominally covers 3 orders of magnitude. Furthermore, as local variations on the order of 10^{−2} m affect the hydrostatic pressure gradient, the solution of η requires precision over at least 5 orders of magnitude – i.e., a stiff numerical solution that can be difficult to converge for either a linear or nonlinear solver. Thus, splitting $\partial \mathit{\eta}/\partial x$ into a downslope body force (S_{R}) and a local residual ($\partial {h}_{\mathrm{a}}/\partial x)$ is effectively removing a largescale gradient from the solution variables, which will generally improve numerical behavior.
Despite the above disadvantages, the η form retains some advantages in creating conservative finitevolume formulations of the SaintVenant equations (Hodges, 2019). Arguably, such methods should be confined to explicit timemarching schemes or localized solutions, in which η covers a smaller range, and the RS method should be preferred for larger systems.
5.3 RS advantages and limitations
The fundamental difference between the SPRNTRS approach and most, if not all, conventional models (including unsteady HECRAS) is that our method algebraically revises the SaintVenant equations to exactly accommodate discontinuous geometry while maintaining a smooth source term, whereas other models typically introduce ad hoc changes (e.g., linearization) to provide stable and faster numerical behavior when discontinuities are likely to cause numerical instabilities. These differences in the governing equations can be expected to cause differences in the solution – especially where nonlinear terms are strong. The present scope is limited in that only a single model was modified (SPRNT/SPRNTRS) and only a single model (HECRAS) was used as an external control. The validity of the underlying algebraic transformation in the RS method has been demonstrated by these tests; however, it remains to be seen how implementing the RS approach in other models – particularly wellbalanced models – might alter residual errors, convergence rates, and computational performance. We are interested in collaborating with other researchers who have access to and familiarity with source code of candidate wellbalanced models.
An important limitation to the present work is that we focus solely on subcritical flow. The Preissmann scheme used in the underlying SPRNT model is known to exhibit instabilities with transcritical flows (Samuels and Skeels, 1990; Sart et al., 2010; Meselhe and Holly, 1997), which can be suppressed with the ad hoc local partial inertia (LPI) scheme of Fread et al. (1996). Our preliminary work (not shown) indicates that the RS approach can stabilize the Preissmann scheme without using LPI, but further work is required to test and validate the RS method for transcritical and supercritical flows.
Overall, the RS method can ensure the Lipschitz smoothness of slope representation in the momentum source term (without smoothing geometry), thus reducing one source of oscillatory or unstable behavior in numerical solutions. However, application of the RS method is not without some limitations. Although the switch from S_{0} to S_{R} is algebraically exact, the application of the RS method requires some method to select the distributed z_{R}(x) and to determine S_{R}(x). Poor selection of z_{R} can theoretically result in nonsmooth S_{R}. In the present work, the profile of z_{R} in the Waller Creek case is generated by the cubic Bspline technique, which is controlled by the number of knots and their spacing. In general, the distance between knots must be longer than the spacing between cross sections so that the generated S_{R} is smooth at the model's discretization scale. It is not clear that a mathematical “optimum” for selection of knots will necessarily exist, but there are likely (unknown) practical limits on knot selection spacing for “adequate” smoothness of z_{R}(x). Our results indicate that approximating cubic Bsplines is adequate for producing smooth z_{R} for the tested geometries and that the solutions are robust for the selection of z_{R} as long as S_{R} is smooth (Yu et al., 2019b). However, it is likely there are limitations to applying the RS method in largescale river network simulations that will make it difficult to use a simple globally applied knot spacing. Such networks might consist of 10^{4} to 10^{5} reaches spanning wide geographical regions with a variety of topology and inconsistent data availability. Some reaches may have welldefined cross sections at close spacing, other reaches might be poorly documented (Hodges, 2013). Thus, it seems likely that a method for automatically generating approximating splines (or some other form of smoothing) would be useful, but such an advance arguably requires a method for quantitatively evaluating the goodness of a particular set of z_{R}(x), which remains an open question. We speculate that simple window filtering techniques may be adequate for river databases such as NHDplus (United States National Hydrography Dataset Plus), but further investigation and examination are needed to better understand the interplay between the smoothing scales and the numerical solution using the RS method for large networks.
This study has implemented RS only in the SPRNT code, as discussed in more detail in Sect. 3. The baseline governing equations for SPRNT are of the form of Eq. (2), the socalled “nonconservative” form, which simply means that the entirety of the hydrostatic pressure gradient is effectively a source term, in contrast to “conservative” equations such as the Cunge–Liggett form (Cunge et al., 1980), in which a portion of the hydrostatic pressure gradient is abstracted to the advection term. Although it remains to be shown in future work, the algebraic transformation implied in Eq. (4) can arguably be applied in the Cunge–Liggett form or any other conservative form of the SVEs. Similarly, the fundamental algebraic transformation to S_{R} and $\partial {h}_{\mathrm{a}}/\partial x$ will be equally valid in any finitevolume method using S_{0} and h_{0}.
The greatest barrier to the adoption of the RS method in an existing model is likely the need to rewrite the geometry functions to accept h_{a}, h_{R}, and z_{R} in place of h_{0} and z_{0}. The difficulties involved in this effort depend on whether or not the model geometry functions are sufficiently isolated from the main solution algorithm. Indeed, we can imagine codes where the geometry functions are essentially dispersed throughout and require extensive effort to alter, debug, and verify.
5.4 The future for RS methods
The RS method as introduced above might be just a starting point. Although the present work focused on the nonconservative form, the concepts presented herein will likely be effective in addressing the wellbalanced problem for conservative forms as reviewed in Kesserwani (2013) and Hodges (2019). Furthermore, the algebra in the RS demonstration above leads to the conjecture that the method could be extended to 2D reference slopes for bathymetry in 2D or 3D models. Undoubtedly there are unknown numerical challenges in extending to higher dimensions – particularly in ensuring a 2D spline function is adequately spaced to ensure smoothness – but there does not appear to be any fundamental conceptual difficulty in such efforts.
The reference slope (RS) method introduces a new form of the SaintVenant equations for 1D river flow. The advantage of the RS method is that it ensures the body force (slope) source term is smooth and cannot destabilize the numerical solution. The RS method introduces the concept of an arbitrary smooth reference elevation, z_{R}(x), with computed reference slopes, S_{R}(x), and associated depths, h_{a}(x). These geometries are algebraically related to the traditional channel thalweg elevation, depth, and bottom slope (${z}_{\mathrm{0}},{h}_{\mathrm{0}},{S}_{\mathrm{0}}$) used in many models. The RS method is implemented in an opensource SaintVenant solver as SPRNTRS. In this study, SPRNTRS was compared to both analytical solutions and the conventional HECRAS model for synthetic test reaches and an urban creek for subcritical flows. The model–model comparisons are within the expected truncation error for both the analytical and synthetic test cases and within acceptable differences for simulating flow through the complex geometry of an urban creek. The slightly larger simulation differences in the urban creek test case are likely due to ad hoc linearization algorithms used in HECRAS that do not appear in SPRNTRS.
As discussed in the Sect. 2, when faced with nonsmooth geometries in a channel reach, prior researchers have resorted to limiting or smoothing discontinuous source terms or employing numerical techniques that mitigate oscillatory/unstable numerical behaviors. In contrast, the new RS method transforms a discontinuous bottom slope source term into a smooth expression without losing either complexity in the geometry or introducing ad hoc smoothing of the geometry, the numerical method, or the solution. An important advantage of the RS method is that it is entirely mechanical, requiring only the selection of control knot spacing for the approximating spline at some length scale larger than the crosssection spacing. That is to say that RS does not require the model designer or user to introduce smoothing thresholds or ad hoc algorithm bounds. As such, we believe the RS method could be particularly valuable as we move from fineresolution reachscale modeling to largescale continental river dynamic simulations (Hodges, 2013) or develop massively parallel stormwater network models for megacities (MoralesHernandez et al., 2020)
The RS method is not specific to SPRNT but can be adapted to any SaintVenant solver that uses a bottom slope (S_{0}) term in the discretization. The mathematical change is conceptually trivial, but the actual effort depends on how crosssection geometry is embedded in the code. The code for both SPRNT and SPRNTRS are available under opensource license at GitHub (Liu, 2014).
As discussed in Sect. 3, the stability of the SPRNTRS simulations for Waller Creek was ensured by merging 36 computational elements in which the crosssection spacing was closer than 10 m. This minimum spacing cutoff was selected as being well below the median spacing of 28.6 m and mean spacing of 37.7 m and proved adequate for ensuring SPRNTRS stability in the tested simulations. Unfortunately, stability of unsteady HECRAS required the further removal of three crosssection elements (shown in Fig. A1) and reducing Manning's n at six additional cross sections (listed in Table A1). Selecting these changes was a matter of art rather than science as we could not identify clear criteria for crosssection removal or Manning's n adjustments for HECRAS – other than that these locations appeared to be where instabilities appeared in unsteady HECRAS model runs. Although SPRNTRS could run without these changes, for consistency in the model comparisons the geometry of the SPRNTRS model was modified to exactly match the adjusted geometry required for the HECRAS model.
A  crosssection area (m^{2}) 
g  gravitational acceleration (ms^{−2}) 
h_{0}  water depth (m) 
h_{a}  associated water depth (m) 
h_{R}  reference height (m) 
n  Gauckler–Manning–Strickler roughness (${\text{m}}^{\mathrm{1}/\mathrm{3}}\text{s}$) 
P_{w}  wetted perimeter (m) 
Q  volumetric flow rate (m^{3}s^{−1}) 
q_{ℓ}  flow rate per unit length through channel sides (m^{2} s^{−1}) 
S_{0}  channel bottom slope 
S_{f}  channel friction slope 
S_{R}  channel reference slope 
S_{SW}  channel sidewall slope 
t  time (s) 
W  channel width (m) 
W_{B}  channel bottom width ((m)) 
x  alongchannel spatial coordinate 
z_{0}  channel bottom elevation (m) 
z_{R}  reference elevation (m) 
α  bottom displacement coefficient 
ϵ  second derivative of water surface elevation (m^{−1}) 
η  water surface elevation (m) 
ρ  normalized difference between results 
ζ  absolute mean normalized difference (AMND) 
The complete code for the reference slope module and SPRNT is available at Github (https://github.com/frankyliu/SPRNT, https://doi.org/10.5281/zenodo.3978737, Liu and Yu, 2020).
All test case files and results have been uploaded to a public repository under the Texas Data Repository (Yu et al., 2019c) and can be accessed through https://doi.org/10.18738/T8/BXJBF5.
CWY and BRH designed and performed the experiments, analyzed the data, and wrote the paper. FL wrote the code for the reference slope method.
The authors declare that they have no conflict of interest.
This research has not been formally reviewed by the U.S. Environmental Protection Agency (EPA). The views expressed in this document are solely those of the authors and do not necessarily reflect those of the agency. The EPA does not endorse any products or commercial services mentioned in this publication.
The authors would like to thank Fernando R. Salas from the U.S. National Water Center for providing assistance on collecting and processing the data for the test cases. The first author, CWY, would also like to thank Alan Plummer Associates, Inc., for providing additional support during the latter stages of this project. Frank Liu would like to acknowledge the support from the Laboratory Directed Research and Development program of Oak Ridge National Laboratory, managed by UTBattelle, LLC, for the U.S. Department of Energy.
This research has been supported by the U.S. National Science Foundation (grant no. CCF1331610) and the U.S. Environmental Protection Agency (cooperative agreement no. 83595001).
This paper was edited by Roberto Greco and reviewed by two anonymous referees.
Aggett, G. R. and Wilson, J. P.: Creating and coupling a highresolution DTM with a 1D hydraulic model in a GIS for scenariobased assessment of avulsion hazard in a gravelbed river, Geomorphology, 113, 21–34, https://doi.org/10.1016/j.geomorph.2009.06.034, 2009. a
Brunner, G. W.: HECRAS, River Analysis System Reference Manual, U.S. Army Corps of Engineers, Hydrologic Engineering Center, Technical Report CPD69, Davis, California, USA, available at: [13:43] Viola Zierenberg https://www.hec.usace.army.mil/software/hecras/documentation/HECRAS 5.0 Reference Manual.pdf (last access: 10 August 2020), 2016a. a, b, c
Brunner, G. W.: HECRAS, River Analysis System User's Manual, Version 5.0, U.S. Army Corps of Engineers, Hydrologic Engineering Center, Technical Report CPD68, Davis, California, USA, available at: https://www.hec.usace.army.mil/software/hecras/documentation/HECRAS 5.0 Users Manual.pdf (last access: 10 August 2020), 2016b. a, b
Burguete, J., GarciaNavarro, P., Murillo, J., and GarciaPalacin, I.: Analysis of the Friction Term in the OneDimensional ShallowWater Model, Journal of Hydraul. Eng.ASCE, 133, 1048–1063, https://doi.org/10.1061/(ASCE)07339429(2007)133:9(1048), 2007. a
Chen, Y. H.: Mathematical Modeling of Water and Sediment Routing in Natural Channels, PhD thesis, Colorado State University, Ft. Collins, CO, 1973. a
Cunge, J. A., Holly, F. M., and Verwey, A.: Practical Aspects of Computational River Hydraulics, Pitman Publishing Ltd, Boston, MA, USA, 1980. a
de Boor, C.: A Practical Guide to Splines, SpringerVerlag, New York Berlin Heidelberg, 2001. a
Decoene, A., Bonaventura, L., Miglio, E., and Saleri, F.: Asymptotic derivation of the sectionaveraged shallow water equations for natural river hydraulics, Math. Mod. Meth. Appl. S., 19, 387–417, https://doi.org/10.1142/S0218202509003474, 2009. a
Di Baldassarre, G.: Floods in a Changing Climate: Inundation Modelling, International Hydrology Series, Cambridge University Press, Cambridge, UK, https://doi.org/10.1017/CBO9781139088411, 2012. a
Fread, D. L., Jin, M., and Lewis, J. M.: An LPI numerical implicit solution for unsteady mixedflow simulation, in: North American Water and Environment Congress, vol. 96, pp. 49–72, 1996. a, b
Gichamo, T. Z., Popescu, I., Jonoski, A., and Solomatine, D.: River crosssection extraction from the ASTER global DEM for flood modeling, Environ. Modell. Softw., 31, 37–46, https://doi.org/10.1016/j.envsoft.2011.12.003, 2012. a
Giustarini, L., Matgen, P., Hostache, R., Montanari, M., Plaza, D., Pauwels, V. R. N., De Lannoy, G. J. M., De Keyser, R., Pfister, L., Hoffmann, L., and Savenije, H. H. G.: Assimilating SARderived water level data into a hydraulic model: a case study, Hydrol. Earth Syst. Sci., 15, 2349–2365, https://doi.org/10.5194/hess1523492011, 2011. a
Greenberg, J. M. and LeRoux, A.Y.: A wellbalanced scheme for the numerical processing of source terms in hyperbolic equations, SIAM Journal on Numerical Analysis, 33, 1–16, https://doi.org/10.1137/0733001, 1996. a
Hicks, F. E. and Peacock, T.: Suitability of HECRAS for flood forecasting, Can. Water Resour. J., 30, 159–174, https://doi.org/10.4296/cwrj3002159, 2005. a
Hodges, B. R.: Challenges in Continental River Dynamics, Environ. Modell. Softw., 50, 16–20, https://doi.org/10.1016/j.envsoft.2013.08.010, 2013. a, b, c
Hodges, B. R.: Conservative finitevolume forms of the SaintVenant equations for hydrology and urban drainage, Hydrol. Earth Syst. Sci., 23, 1281–1304, https://doi.org/10.5194/hess2312812019, 2019. a, b, c
Horritt, M. S. and Bates, P. D.: Evaluation of 1D and 2D numerical models for predicting river flood inundation, J. Hydrol., 268, 87–99, https://doi.org/10.1016/S00221694(02)00121X, 2002. a
Iserles, A.: A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, Cambridge, UK, 1996. a
Kesserwani, G.: Topography discretization techniques for Godunovtype shallow water numerical models: a comparative study, J. Hydraul. Res., 51, 351–367, https://doi.org/10.1080/00221686.2013.796574, 2013. a, b, c, d
Kuiry, S. N., Sen, D., and Bates, P. D.: Coupled 1DQuasi2D Flood Inundation Model with Unstructured Grids, J. Hydraul. Eng.ASCE, 136, 493–506, 2010. a
Liang, Q. and Marche, F.: Numerical resolution of wellbalanced shallow water equations with complex source terms, Adv. Water Resour., 32, 873–884, https://doi.org/10.1016/j.advwatres.2009.02.010, 2009. a, b
Liggett, J. A.: Numerical method of solution of the unsteady flow equations, Unsteady Flow in Open Channels, 1, 89–182, 1975. a
Liu, F.: SPRNT: A River Dynamics Simulator, available at: https://github.com/frankyliu/SPRNT, last access; 4 September, 2018, 2014. a, b, c
Liu, F. and Hodges, B. R.: Applying microprocessor analysis methods to river network modeling, Environ. Modell. Softw., 52, 234–252, https://doi.org/10.1016/j.envsoft.2013.09.013, 2014. a, b, c
Liu, F. and Yu, C.W.: frankyliu/SPRNT: minor bug fixes (Version v1.3.8), Zenodo, https://doi.org/10.5281/zenodo.3978737, 2020. a
MacDonald, I., Baines, M., Nichols, N., and Samuels, P.: Comparison of some steady state SaintVenant solvers for some test problems with analytic solutions, Numerical Analysis Report 2/95, Department of Mathematics, University of Reading, available at: http://www.reading.ac.uk/web/files/maths/0295.pdf (last access: 10 August 2020), 1995. a, b, c, d, e, f, g, h, i, j, k, l
MartinezAranda, S., Murrillo, J., and GarciaNavarro, P.: A 1D numerical model for the simulation of unsteady and highly erosive flows in rivers, Comput. Fluids, 181, 8–34, https://doi.org/10.1016/j.compfluid.2019.01.011, 2019. a
Mejia, A. I. and Reed, S. M.: Evaluating the effects of parameterized cross section shapes and simplified routing with a coupled distributed hydrologic and hydraulic model, J. Hydrol., 409, 512–524, https://doi.org/10.1016/j.jhydrol.2011.08.050, 2011. a
Meselhe, E. and Holly Jr., F.: Invalidity of Preissmann scheme for transcritical flow, J. Hydraul. Eng., 123, 652–655, https://doi.org/10.1061/(ASCE)07339429(1997)123:7(652), 1997. a
MoralesHernández, M., Sharif, M. B., Gangrade, S., Dullo, T. T. , Kao, S.C., Kalyanapu, A., Ghafoor, S. K., Evans, K. J., MadadiKandjani, E., and Hodges, B. R.: Highperformance computing in water resources hydrodynamics. Journal of Hydroinformatics, jh2020163, https://doi.org/10.2166/hydro.2020.163, 2020. a
Nujic, M.: Efficient implementation of nonoscillatory schemes for the computation of freesurface flows, J. Hydraul. Res., 33, 101–111, https://doi.org/10.1080/00221689509498687, 1995. a
Preissmann, A.: Progagation des intumescences dans les canaux et rivières, in: Proceedings of First Congress of French Association for Computation, Grenoble, France, 433–442, 1961. a
Samuels, P. G. and Skeels, C. P.: Stability limits for Preissmann's scheme, J. Hydraul. Eng., 116, 997–1012, https://doi.org/10.1061/(ASCE)07339429(1990)116:8(997), 1990. a
Sanders, B. F.: Highresolution and nonoscillatory solution of the St. Venant equations in nonrectangular and nonprismatic channels, J. Hydraul. Res., 39, 321–330, https://doi.org/10.1080/00221680109499835, 2001. a, b
Sart, C., Baume, J.P., Malaterre, P.O., and Guinot, V.: Adaptation of Preissmann's scheme for transcritical open channel flows, J. Hydraul. Res., 48, 428–440, https://doi.org/10.1080/00221686.2010.491648, 2010. a
Schippa, L. and Pavan, S.: Analytical treatment of source terms for complex channel geometry, J. Hydraul. Res., 46, 753–763, 2008. a
Sharkey, J. K.: Investigating Instabilities with HECRAS Unsteady Flow Modeling for Regulated Rivers at Low Flow Stages, MS thesis, University of Tennessee, Knoxville, Tennessee, USA, available at: https://trace.tennessee.edu/utk_gradthes/3183 (last access: 10 August 2020), 2014. a
Singh, J., Altinakar, M. S., and Ding, Y.: Numerical modeling of rainfallgenerated overland flow using nonlinear shallowwater equations, J. Hydrol. Eng., 20, 04014089, https://doi.org/10.1061/(ASCE)HE.19435584.0001124, 2015. a
Song, L., Zhou, J., Li, Q., Yang, X., and Zhang, Y.: An unstructured finite volume model for dambreak floods with wet/dry fronts over complex topography, Int. J. Numer. Meth. Fl., 67, 960–980, 2011. a
Tayfur, G., Kavvas, M. L., Govindaraju, R. S., and Storm, D. E.: Applicability of St. Venant equations for twodimensional overland flows over rough infiltrating surfaces, J. Hydraul. Eng., 119, 51–63, https://doi.org/10.1061/(ASCE)07339429(1993)119:1(51), 1993. a
Tseng, M.H.: Improved treatment of source terms in TVD scheme for shallow water equations, Adv. Water Resour., 27, 617–629, https://doi.org/10.1016/j.advwatres.2004.02.023, 2004. a
Wang, W., Yang, X., and Yao, T.: Evaluation of ASTER GDEM and SRTM and their suitability in hydraulic modelling of a glacial lake outburst flood in southeast Tibet, Hydrol. Process., 26, 213–225, https://doi.org/10.1002/hyp.8127, 2012. a
Xia, X. and Liang, Q.: A new efficient implicit scheme for discretising the stiff friction terms in the shallow water equations, Adv. Water Resour., 117, 87–97, 2018. a, b
Yu, C.W., Liu, F., and Hodges, B. R.: Consistent initial conditions for the SaintVenant equations in river network modeling, Hydrol. Earth Syst. Sci., 21, 4959–4972, https://doi.org/10.5194/hess2149592017, 2017. a
Yu, C.W., Hodges, B. R., and Liu, F.: Evaluation of the Performance of the Reference Slope (RS) Method in Simulation Program for River Network (SPRNT), Texas Data Repository Dataverse, https://doi.org/10.18738/T8/GRCB8B, 2019a. a, b
Yu, C.W., Hodges, B. R., and Liu, F.: Evaluation of the Performance of the Reference Slope (RS) Method in Simulation Program for River Network (SPRNT), Tech. rep., University of Texas at Austin, https://doi.org/10.18738/T8/WMQIPS, 2019b. a, b
Yu, C.W., Hodges, B. R., and Liu, F.: Test case data for Reference Slope study with HECRAS and SPRNTRS, Texas Data Repository Dataverse, V1, https://doi.org/10.18738/T8/BXJBF5, 2019c. a
Zhou, J. G., Causon, D. M., Mingham, C. G., and Ingram, D. M.: The surface gradient method for the treatment of source terms in the shallowwater equations, J. Comput. Phys., 168, 1–25, https://doi.org/10.1006/jcph.2000.6670, 2001. a
 Abstract
 Introduction
 Background
 Methods
 Results
 Discussion
 Conclusions
 Appendix A: Geometry adjustments for stable unsteady HECRAS simulations
 Appendix B: Notation
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Background
 Methods
 Results
 Discussion
 Conclusions
 Appendix A: Geometry adjustments for stable unsteady HECRAS simulations
 Appendix B: Notation
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References