A new form of the Saint-Venant equations for variable topography

Abstract

The solution stability of river models using the one-dimensional (1D) Saint-Venant 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 Saint-Venant 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 cross-section 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 real-world geometry. The validity of the mathematical concept is demonstrated with the open-source 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 (HEC-RAS) model. The new method reduces numerical oscillations and instabilities without requiring ad hoc smoothing algorithms.

1 Introduction

The Saint-Venant equations (SVEs) for one-dimensional (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_{\mathrm{R}}$) and its associated depth ($h_{\mathrm{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_{\mathrm{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_{\mathrm{R}}$, to distinguish it from the traditional thalweg bottom slope, $S_{\mathrm{0}}$.

The two common differential forms of the SVEs are

1$${\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)=-gA{S}_{\mathrm{f}}-gA{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}$$ and 2$${\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)=-gA{S}_{\mathrm{f}}+gA{S}_{\mathrm{0}}-gA{\displaystyle \frac{\partial h_{\mathrm{0}}}{\partial x}},$$ where $Q$ is the flow rate, $A$ is the cross-section area, $\mathit{\eta }$ is the water surface elevation, $h_{\mathrm{0}}$ is the thalweg depth, $S_{\mathrm{0}}$ is the thalweg bottom slope, and $S_{\mathrm{f}}$ is the friction slope that represents the local energy gradient. Equation () can be envisioned as using the Piezometric head gradient to force the flow, as shown in Fig. a. In contrast, Eq. () 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. b. 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.

Figure 1

The $\mathit{\eta }$ form of the SVEs (a) has a driving Piezometric head gradient which is equivalent (b) to the sum of the hydrostatic head gradient and a body force aligned with $S_{\mathrm{0}}$. The effect of varying geometry is handled in $A$ in both forms.

[Figure omitted. See PDF]

The two SVE forms of Eqs. () and () are algebraically identical using the identity 3$${\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}\equiv {\displaystyle \frac{\partial h_{\mathrm{0}}}{\partial x}}-S_{\mathrm{0}}.$$ Heuristically, we can propose a more general identity of 4$${\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}\equiv {\displaystyle \frac{\partial h_{\mathrm{a}}}{\partial x}}-S_{\mathrm{R}},$$ where $S_{\mathrm{R}}$ is an arbitrary reference slope and $h_{\mathrm{a}}$ is an associated depth that will be defined in Sect.  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_{\mathrm{R}}$. Applying Eq. () to () provides the following: 5$${\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)=-gA{S}_{\mathrm{f}}+gA{S}_{\mathrm{R}}-gA{\displaystyle \frac{\partial h_{\mathrm{a}}}{\partial x}},$$ where the terms $gA{S}_{\mathrm{R}}$ and $gA\frac{\partial h_{\mathrm{a}}}{\partial x}$ are algebraically equivalent to those in Eqs. () and (). Clearly, if we let $S_{\mathrm{R}}$ equal $S_{\mathrm{0}}$, then $h_{\mathrm{a}}$ equals $h_{\mathrm{0}}$ and we recover Eq. (). Furthermore, if we let $S_{\mathrm{R}}$ equal 0, then $h_{\mathrm{a}}$ equals $\mathit{\eta }$ and we recover Eq. (). The equations are identical with these substitutions, so it follows that using a reference slope of zero ($S_{\mathrm{R}}=\mathrm{0}$) must exactly represent the same topographic variability as using a reference slope that mimics the topographic slope ($S_{\mathrm{R}}=S_{\mathrm{0}}$) as long as $h_{\mathrm{a}}$ is correctly defined consistent with Eq. (). That is to say that from simple algebra the use of the real $S_{\mathrm{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 cross-section flow areas are correctly computed.

From the arguments above, the effects of varying bottom topography are captured by $S_{\mathrm{R}}=\mathrm{0}$ and $h_{\mathrm{a}}=\mathit{\eta }$, which implies we are also free to introduce any other (preferably smooth) $S_{\mathrm{R}}$ into Eq. () without altering the underlying representation of variable topography. An example is illustrated in Fig. . As the splitting defined in Eq. () makes Eq. () algebraically identical to Eqs. () and (), the introduction of a smooth $S_{\mathrm{R}}$ does not reflect “smoothing” of the topography. It merely reflects a decision on whether effects of non-smoothness will reside solely in solution variables $A$ and $h_{\mathrm{a}}$ or will also be forced as a non-smooth source term in $S_{\mathrm{0}}$.

Figure 2

Comparison of Piezometric forcing terms for varying topography. (a) Eq. () with the arbitrary $S_{\mathrm{R}}$ line provided for reference. (b) The equivalent split form of forcing with Eq. () using the identity of Eq. (). The physical bottom topography, shown only in (a) for clarity, only plays a role through the cross-section area ($A$) which feeds back into the solution of variability in $h_{\mathrm{a}}$ in both forms of the equation.

[Figure omitted. See PDF]

We would like to use an a priori smooth $S_{\mathrm{R}}$ in a computational model rather than the actual thalweg $S_{\mathrm{0}}$ because of what happens to $S_{\mathrm{0}}\left(x\right)$ and $\partial h_{\mathrm{0}}\left(x\right)/\partial x$ for topography varying sharply over short distances, as illustrated in Fig. . From a physics perspective, using $S_{\mathrm{0}}$ to split the Piezometric head is an intuitive way to describe the local interplay of pressure with the bottom slope. Furthermore, $S_{\mathrm{0}}$ has the advantage of being readily reduced to a kinematic wave equation where $S_{\mathrm{f}}=S_{\mathrm{0}}$, which has some advantage in multipurpose codes. However, from a numerical modeling perspective, using $S_{\mathrm{0}}$ has a significant limitation based on its smoothness. If the water surface is smooth, then non-smooth $S_{\mathrm{0}}\left(x\right)$ requires the numerical solver to produce a compensating non-smooth $h_{\mathrm{0}}\left(x\right)$, i.e., requiring a well-balanced scheme (see Sect. ). If we can discard our (wrong) intuition that the $S_{\mathrm{0}}$ form must somehow better represent sharply variable topography – i.e., recognizing the algebraic equivalence of Eq. () with Eqs. () and () – it follows that the splitting of the Piezometric head to include a body force that is everywhere exactly aligned with a sharply varying $S_{\mathrm{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_{\mathrm{R}}$ to be smooth, we can ensure the $h_{\mathrm{a}}$ solution is also smooth for a smoothly varying free surface.

Figure 3

Comparison of Piezometric forcing terms with sharply varying topography for Eqs. () and (). (a) The Piezometric head gradient for Eq. (). (b) The equivalent split form of forcing from Eq. () using the identity of Eq. (). The variability of the bottom of (a) will be reflected in the cross-section area, $A\left(x\right)$, that feeds back into the solution of $\mathit{\eta }\left(x\right)$. In contrast, (b) attempts to directly represent topographic variability as both a driving term in $S_{\mathrm{0}}\left(x\right)$ and a solution response term, $h_{\mathrm{0}}\left(x\right)$

[Figure omitted. See PDF]

The use of $S_{\mathrm{R}}$ rather than $S_{\mathrm{0}}$ in the governing equations can perhaps be better understood if we think of the slope in Eq. () 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.  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. .

In this paper, we examine the effect of variable cross-section 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_{\mathrm{0}}$ is essentially unchanged. We further demonstrate that bottom slope discontinuities are a cause of problems in finite-difference forms of 1D Saint-Venant 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_{\mathrm{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_{\mathrm{R}}\left(x\right)$, which results in smooth reference slopes, $S_{\mathrm{R}}\left(x\right)$, and the redefinition of the thalweg depth ($h_{\mathrm{0}}$) as a depth associated ($h_{\mathrm{a}}$) with the reference elevation. Note that $h_{\mathrm{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 Saint-Venant model as it is entirely independent of the solution algorithm; however, implementation does require rewriting code for the relationships between cross-section 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_{\mathrm{0}}$) as the depth associated ($h_{\mathrm{a}}$) with the reference elevation, $z_{\mathrm{R}}$.

2 Background

One-dimensional (1D) hydrodynamic models using the Saint-Venant equations (SVEs) are widely employed for studying both natural streams and man-made channels e.g.,. It is widely recognized that numerical solutions of the SVEs are prone to spurious oscillations in the free-surface elevation unless particular care is taken in the numerical formulation and/or the problem definition e.g.,. Numerous techniques and special numerical schemes have been previously designed to overcome unwanted numerical oscillations caused by discontinuous geometries and boundary conditions e.g.,. These approaches typically rely on the concept of a well-balanced discrete form as discussed in a comprehensive review by and further elaborated by .

Although well-balanced schemes are relatively robust in handling discontinuous boundary conditions, they have not been extensively applied in water resource models to simulate the regional-to-continental-scale river networks or stormwater systems for megacities. The rapidly varying and discontinuous $S_{\mathrm{0}}$ in natural systems can significantly increase the difficulty and computational burden of obtaining a well-balanced method . Hence, when a large-scale open-channel 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 . 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 non-smooth source terms; i.e., if we consider an advection equation of the form

6$${\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)=\mathit{\sigma },$$ where $\mathit{\sigma }$ 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.,. That numerical instabilities are often caused by non-smooth 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 , positivity-preserving methods for coupled models , and implicit schemes that address stiffness of the nonlinear friction term . 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_{\mathrm{0}}$, whose non-smoothness 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 well-balanced schemes (see reviews noted above) seek to compensate for non-smoothness of all parts of the source term in the structure of the numerical discretization. Arguably, if the slope term is guaranteed smooth, then a well-balanced scheme should be simpler to create. In general, when the thalweg bottom slope ($S_{\mathrm{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_{\mathrm{0}}\left(x\right)$ 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_{\mathrm{0}}$ will tend to have significant variability. It follows that the use of $S_{\mathrm{0}}$ has the undesirable property that smaller $\mathrm{\Delta }x$ (i.e., resolving a river with more detailed survey data) will increase the non-smoothness 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_{\mathrm{0}}$ smoothness, when it occurs in a model of a natural river channel, is typically the result of relatively long separations ($\mathrm{\Delta }x)$ between cross-section surveys that ensures that the discrete $d^{\mathrm{2}}{z}_{\mathrm{0}}/\mathrm{d}{x}^{\mathrm{2}}$ is small. Thus, removing cross sections can be an effective mitigation technique because it increases $\mathrm{\Delta }x$ and effectively smooths $S_{\mathrm{0}}$. In general, models discretized with higher-resolution river surveys (smaller $\mathrm{\Delta }x$) will have greater non-smoothness in $S_{\mathrm{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_{\mathrm{0}}$ can be understood by considering the identity in Eq. () for a channel with subcritical flow in which the free-surface curvature is expected to be negligible, i.e., ${\partial }^{\mathrm{2}}\mathit{\eta }/\partial x^{\mathrm{2}}\sim \mathit{ϵ}$, where $\mathit{ϵ}\ll \text{d}{S}_{\mathrm{0}}/\text{d}x,{\partial }^{\mathrm{2}}\mathit{\eta }/{\partial }^{\mathrm{2}}x$. Taking the along-channel gradient of Eq. () implies that 7$${\displaystyle \frac{{\partial }^{\mathrm{2}}{h}_{\mathrm{0}}}{\partial x^{\mathrm{2}}}}={\displaystyle \frac{\text{d}{S}_{\mathrm{0}}}{\text{d}x}}-\mathit{ϵ}.$$ Thus, forcing with a non-smooth $\text{d}{S}_{\mathrm{0}}\left(x\right)/\text{d}x$ will require non-negligible curvature of the response variable, $h_{\mathrm{0}}\left(x\right)$, whose gradient is also a forcing function of the nonlinear equation. Feedback can easily build and cause successive overshoot/undershoot effects, producing oscillations and non-convergence in a nonlinear solver. In contrast, Eq. () can be invoked with $\text{d}{S}_{\mathrm{R}}/\text{d}x$ 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 non-smooth. As a practical matter, any $S_{\mathrm{0}}$ with a discontinuous discrete first derivative (i.e., discontinuities in the second derivative of the thalweg elevation, ${\text{d}}^{\mathrm{2}}{z}_{\mathrm{0}}/\text{d}{x}^{\mathrm{2}}$) will be Lipschitz discontinuous and should not be directly discretized in an SVE solution with Eq. (). Although approximate numerical solutions of equations with non-smooth $S_{\mathrm{0}}\left(x\right)$ can sometimes be attained for models with sufficient damping, such solutions are questionable as they do not have rigorous mathematical foundations.

Arguably, non-smoothness in $S_{\mathrm{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 non-smoothness – e.g., the well-balanced 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. () is the extreme example of the last approach – replacing $S_{\mathrm{0}}$ and $\partial h_{\mathrm{0}}/\partial x$ with $\partial \mathit{\eta }/\partial x$ ensures that $S_{\mathrm{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. ) 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 free-surface elevation despite non-smooth geometry. Nevertheless, common SVE solvers can exhibit oscillatory, non-convergent 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 non-smooth jumps) necessarily requires more analyses than is practical in a single paper.

3 Methods

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. (), which is coupled to the solution of continuity,

8$${\displaystyle \frac{\partial A}{\partial t}}+{\displaystyle \frac{\partial Q}{\partial x}}=q_{\mathrm{\ell }},$$ where $q_{\mathrm{\ell }}$ is a lateral inflow per unit length. Note that significantly non-smooth $q_{\mathrm{\ell }}(x,t)$ can provide another source of numerical oscillations and instability . 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_{\mathrm{\ell }}=\mathrm{0}$ everywhere.

In the SPRNT momentum equation, the friction slope is represented using the Chezy–Manning form as 9$$S_{\mathrm{f}}={\displaystyle \frac{n^{\mathrm{2}}{P}_{\mathrm{w}}^{\mathrm{4}/\mathrm{3}}}{A^{\mathrm{10}/\mathrm{3}}}}{Q}^{\mathrm{2}},$$ where $P_{\mathrm{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.,, the Chezy–Manning form remains popular due to its simplicity.

The baseline model uses the thalweg elevation ($z_{\mathrm{0}}$), the thalweg depth ($h_{\mathrm{0}}$), and the thalweg bottom slope ($S_{\mathrm{0}}$) as 10$$h_{\mathrm{0}}\equiv \mathit{\eta }-z_{\mathrm{0}}$$ and 11$$S_{\mathrm{0}}\equiv -{\displaystyle \frac{\partial z_{\mathrm{0}}}{\partial x}}.$$ SPRNT is an open-source, 1D hydrodynamic solver using the fully implicit Preissmann numerical scheme with the Newton–Raphson iteration and computational acceleration techniques developed from very large-scale integration (VLSI) semiconductor design. Details on the baseline SPRNT model and its application to large river networks are provided in and .

3.2 Reference slope method

We introduce a new reference slope (RS) method through a transformation and redefinition of geometry in the Saint-Venant equations, as discussed in Sect. . In place of the conventional $h_{\mathrm{0}}$ and $z_{\mathrm{0}}$, we define a reference elevation ($z_{\mathrm{R}}$) and its associated depth ($h_{\mathrm{a}}$) as shown in Fig. . These provide a relationship with the free-surface elevation ($\mathit{\eta }$) defined as

12$$h_{\mathrm{a}}\equiv \mathit{\eta }-z_{\mathrm{R}}.$$ Note that $z_{\mathrm{R}}$ is arbitrary, so $h_{\mathrm{a}}$ may be either greater than or less than the thalweg depth, $h_{\mathrm{0}}$, at a given location. As shown in Fig. , it is convenient to define the reference height, $h_{\mathrm{R}}$, in relation to $z_{\mathrm{R}}$ and the true bottom elevation, $z_{\mathrm{0}}$, by 13$$h_{\mathrm{R}}\equiv z_{\mathrm{0}}-z_{\mathrm{R}}.$$ Thus, the conventional $h_{\mathrm{0}}$ and $z_{\mathrm{0}}$ are recovered with 14$$h_{\mathrm{0}}=h_{\mathrm{a}}-h_{\mathrm{R}}$$ and 15$$z_{\mathrm{0}}=z_{\mathrm{R}}+h_{\mathrm{R}}.$$ We define the reference slope ($S_{\mathrm{R}}$) as the downstream slope of $z_{\mathrm{R}}$: 16$$S_{\mathrm{R}}\equiv -{\displaystyle \frac{\partial z_{\mathrm{R}}}{\partial x}}.$$ Using Eqs. () and () in Eq. () provides Eq. (), which is more conveniently written as 17$${\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)=-gA{\displaystyle \frac{\partial h_{\mathrm{a}}}{\partial x}}+gA\left(S_{\mathrm{R}}-S_{\mathrm{f}}\right).$$ The above is identical to Eq. () with the simple substitution of $h_{\mathrm{a}}$ and $S_{\mathrm{R}}$ for $h_{\mathrm{0}}$ and $S_{\mathrm{0}}$. In this formulation, the definition of $z_{\mathrm{R}}\left(x\right)$ is arbitrary, so we can a priori require a definition such that $S_{\mathrm{R}}\left(x\right)$ is smooth. A trivial choice that is guaranteed smooth is $z_{\mathrm{R}}\left(x\right)=\text{constant}$, which returns $S_{\mathrm{R}}\left(x\right)=\mathrm{0}$ and the Saint-Venant equations in the form of Eq. (). 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 $\mathit{\eta }$, e.g., the elevation change in a river from its mountain source to a coastal plain . Using the conventional $S_{\mathrm{0}}$ in Eq. () reduces this problem as the range of $h_{\mathrm{0}}$ is inherently confined to the local water depths rather than the underlying topography. In the RS method, the range of $h_{\mathrm{a}}$ is tied to the range of water depths and the selection of $z_{\mathrm{R}}$; thus, for present purposes, we are interested in nontrivial definitions of $z_{\mathrm{R}}$ that (i) are close to $z_{\mathrm{0}}$ to maintain a small range of $h_{\mathrm{a}}$ values and (ii) provide smooth $S_{\mathrm{R}}\left(x\right)$. An arbitrary $z_{\mathrm{R}}$ that is far from $z_{\mathrm{0}}$ or non-smooth is of little interest as it holds no theoretical or practical advantage over the Eq. () approach implied by $S_{\mathrm{R}}\left(x\right)=\mathrm{0}$.

If $S_{\mathrm{R}}\left(x\right)$ is required to be smooth, then the source term of the equation can be guaranteed smooth as long as $S_{\mathrm{f}}\left(x\right)$ is smooth, which is typically true as long as the solution, $Q\left(x\right)$, is smooth. Note that in extreme cases of geometric discontinuity the combined values of $n$, $P_{\mathrm{w}}$, and $A$ in Eq. () can cause a non-Lipschitz 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_{\mathrm{w}}$ as they are coefficients of the solution variables $\left\{A,Q\right\}$. More subtle problems might arise due to discontinuities developed in the $Q^{\mathrm{2}}{A}^{-\mathrm{10}/\mathrm{3}}$ ratio in Eq. (); countering incipient instabilities from this term requires other numerical strategies e.g.,.

Figure 4

Relationships of $h_{\mathrm{a}}$, $h_{\mathrm{R}}$, $z_{\mathrm{R}}$, and $z_{\mathrm{0}}$ for an arbitrary cross section. Note that $z_{\mathrm{R}}>z_{\mathrm{0}}$ is also allowed, which results in $h_{\mathrm{a}}<h_{\mathrm{0}}$ and a negative value for $h_{\mathrm{R}}$. Furthermore, if $z_{\mathrm{R}}$ is greater than $z_{\mathrm{0}}+h_{\mathrm{0}}$, then $h_{\mathrm{a}}$ is negative to retain algebraic consistency. Modified from and used with permission.

[Figure omitted. See PDF]

A critical change required by the introduction of $h_{\mathrm{a}}$ is that the conventional geometric auxiliary relationships of $A=f\left(h_{\mathrm{0}}\right)$ and $P_{\mathrm{w}}=f\left(h_{\mathrm{0}}\right)$ must be transformed into $A=f\left(h_{\mathrm{a}}\right)$ and $P_{\mathrm{w}}=f\left(h_{\mathrm{a}}\right)$. That is, once we change the depth in our equation from $h_{\mathrm{0}}$ to $h_{\mathrm{a}}$, we must re-index the geometry. In general, for known functions $A\left(h_{\mathrm{0}}\right)$ and $P_{\mathrm{w}}\left(h_{\mathrm{0}}\right)$ this is a trivial transformation: 18$$A=A(h_{\mathrm{a}}-h_{\mathrm{R}})$$ and 19$$P_{\mathrm{w}}=P_{\mathrm{w}}(h_{\mathrm{a}}-h_{\mathrm{R}}).$$ However, the implementation in an existing code is not necessarily as simple as the above equations suggest. For example, in Fig. , the condition $A=\mathrm{0}$ occurs when $h_{\mathrm{a}}=h_{\mathrm{R}}$, i.e., a nonzero value as compared to $h_{\mathrm{0}}=\mathrm{0}$ with conventional geometry. Unfortunately, model developers typically have ad hoc wetting/drying treatments that are introduced as $h_{\mathrm{0}}\to \mathrm{0}$ or for $h_{\mathrm{0}}<\mathrm{0}$. Such treatments need to be modified to be deployed as $h_{\mathrm{a}}\to h_{\mathrm{R}}$, which introduces the added complication that $h_{\mathrm{R}}$ is negative when $z_{\mathrm{R}}>z_{\mathrm{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 well-balanced numerical discretization scheme for the conventional $S_{\mathrm{0}}$ representation of geometry e.g.,.

The modification of the SPRNT code to implement the above RS method will be known as SPRNT-RS. The SPRNT and SPRNT-RS source codes are available in an open-source repository . Note that the solution algorithm for SPRNT-RS is identical to that of SPRNT; the only code changes are for the new geometry definitions for $h_{\mathrm{a}}$ and $S_{\mathrm{R}}$ that replace $h_{\mathrm{0}}$ and $S_{\mathrm{0}}$ in the original algorithm. This simple geometry replacement strategy is effective because Eq. () is identical to Eq. () except for the change in nomenclature to $h_{\mathrm{a}}$ and $S_{\mathrm{R}}$.

3.3

Generating a smooth $S_{\mathrm{R}}\left(x\right)$

The $z_{\mathrm{R}}\left(x\right)$ and hence $S_{\mathrm{R}}\left(x\right)$ are arbitrary choices in the RS method but should be generated for the smoothness of $S_{\mathrm{R}}\left(x\right)$ 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_{\mathrm{R}}$ or produced by a different order of splines. For our urban creek test case, the $z_{\mathrm{R}}$ is generated with an approximating cubic B-spline based on thalweg, $z_{\mathrm{0}}\left(x\right)$, elevations. There are a variety of possible ways to generate smooth $S_{\mathrm{R}}\left(x\right)$, but applying approximating cubic B-splines to the $z_{\mathrm{0}}\left(x\right)$ is convenient because the slope is guaranteed to be locally smooth as long as the knot spacing in the B-spline 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 cross-section spacing but non-smooth at the model's discretization scale. That is to say that the exact cubic spline fitting of $z_{\mathrm{0}}\left(x\right)$ does not reduce discontinuities at the discretization scale – only an approximate fitting associated with coarser scales than the cross-section spacing will be effective. It follows that there is some (limited) choice in the selection of the subset of $z_{\mathrm{0}}\left(x\right)$ used as the spline knots, with different sets producing slightly different $\left\{z_{\mathrm{R}}\left(x\right),h_{\mathrm{R}}\left(x\right)\right\}$ 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_{\mathrm{R}}\left(x\right)$ are sufficiently smooth. Implications of the method chosen for generating $z_{\mathrm{R}}\left(x\right)$ are discussed in Sect. . Further details and test cases are provided in .

3.4 HEC-RAS 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 HEC-RAS. showed SPRNT simulations agreed with HEC-RAS with $\le \mathrm{3}\mathit{\%}$ difference in water depth solution when using both prismatic cross sections and nonuniform channels. Thus, HEC-RAS provides a reasonable model for testing and validating SPRINT-RS. We would have preferred to use a single model with and without the RS method for such model–model comparisons; however, HEC-RAS is a closed-source proprietary model, so we could not directly implement and test the RS method in that code. Conversely, as expected by the discussions in Sect. , the baseline SPRNT model is oscillatory and non-convergent on the highly discontinuous geometry of our test cases due to its use of the $S_{\mathrm{0}}$ approach, so it cannot be directly used for before and after comparisons of RS. Thus, simulations using SPRNT-RS are compared to HEC-RAS simulations for validation and insight.

HEC-RAS provides a convenient validation model for three reasons. Firstly, it is a widely accepted engineering model for river-reach simulations, e.g.,. Secondly, it has been used as a validation model in numerous prior studies e.g.,. Finally, unsteady HEC-RAS employs $\partial \mathit{\eta }/\partial x$ as the piezometric gradient rather than using $\partial h_{\mathrm{0}}/\partial x$ and $S_{\mathrm{0}}$, which is one of the reasons it is relatively robust for non-smooth 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 with Lipschitz continuous geometry, various prismatic shapes, and different formulations of $S_{\mathrm{R}}$, (ii) seven synthetic test cases using Lipschitz discontinuous geometry, and (iii) an urban creek with complex cross-section geometry derived from physical surveys that include discontinuities an-order-of-magnitude 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 are used to show that SPRNT-RS reproduces the correct water surface elevation regardless of the selection of $S_{\mathrm{R}}$. These test cases are representative of the more comprehensive analysis provided in . The configuration details for each case are provided in Table , for which we adopt the nomenclature of 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. . To illustrate the adaptability of the RS method, the UR3 and UT2 cases use splines that produce $z_{\mathrm{R}}$ very close to (but not identical to) the actual bed, whereas the other cases use uniform $S_{\mathrm{R}}$ or splines with greater differences. The uniform $S_{\mathrm{R}}$ in cases UR1, UT1, and VR1 is set to the average slope in each domain. With reference to Fig. , the differences between the channel bottom ($z_{\mathrm{0}}$) and the reference bottom ($z_{\mathrm{R}}$) shown in Fig.  imply channel bottom offsets ($h_{\mathrm{R}}$) of varying complexity for the RS method, as shown in Fig. . The VR1 and VR2 cases also provide smooth changes in the channel width, which are shown in Fig. . The node spacing for all these tests is a uniform $\mathrm{\Delta }x=\mathrm{10}$ m. The boundary conditions follow .

Table 1

Configuration and geometric data for analytical test cases derived from . $W_{\mathrm{B}}$ and $S_{\mathrm{SW}}$ represent bottom width and sidewall slope, respectively.

Case name	Cross-section shape type	Cross-section shape detail	$n$	$Q$	$S_{\mathrm{R}}$
UR1	Uniform rectangular	$W_{\mathrm{B}}=\mathrm{10}$ m	0.03	20 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$	Constant $z_{\mathrm{R}}$
UR3	Uniform rectangular				First-order spline
UT1	Uniform trapezoidal	$W_{\mathrm{B}}=\mathrm{9}$ m, $S_{\mathrm{SW}}=\mathrm{2}$			Constant $z_{\mathrm{R}}$
UT2	Uniform trapezoidal	$W_{\mathrm{B}}=\mathrm{10}$ m, $S_{\mathrm{SW}}=\mathrm{2}$			Second-order spline
VR1	Varying rectangular	Varying $W_{\mathrm{B}}$			Constant $z_{\mathrm{R}}$
VR2	Varying rectangular				Third-order spline

Figure 5

Channel bed elevation ($z_{\mathrm{0}}$) and reference bed elevation ($z_{\mathrm{R}}$) for six test cases from .

[Figure omitted. See PDF]

Figure 6

Reference bed offset, $h_{\mathrm{R}}$, for $z_{\mathrm{R}}$ and $z_{\mathrm{0}}$ of the test cases in Fig. .

[Figure omitted. See PDF]

Figure 7

Plan view of channel edges for cases VR1 (a) and VR2 (b) of Fig. .

[Figure omitted. See PDF]

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 non-smooth bathymetry in which most models using $h_{\mathrm{0}}$ and $S_{\mathrm{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 HEC-RAS model for comparison. The simulations use time-invariant boundary conditions with geometry defined by trapezoidal cross sections of uniform side slope, as detailed in Table . The channel bed offset ($h_{\mathrm{R}}$) and thalweg slope ($S_{\mathrm{0}}$) are illustrated in Fig. . 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${}^{\circ }$ sidewall slopes. Bottom roughness is fixed by a Manning's $n$ value of 0.04 for all segments.

Table 2

Configuration of synthetic channel reach test cases: $\mathit{\alpha }$ is used in Eq. () for random perturbation of the baseline Case 1 geometry; $h_{\mathrm{R}}$ is the bed offset based on Fig.  with parentheses indicating upper and lower limits of randomized geometry values over the nonuniform test reach; $S_{\mathrm{0}}$ is the range of the thalweg slope for $z_{\mathrm{0}}$ from Eq. (); and $S_{\mathrm{R}}$ is the selected uniform reference slope.

Case	$\mathit{\alpha }$	$h_{\mathrm{R}}$	$S_{\mathrm{0}}$	$S_{\mathrm{R}}$
Case 1 (baseline)	–	–	0.008	0.008
Case 2	0.01	( $-\mathrm{0.0018}$, $\mathrm{0.0012}$)	(0.0079, 0.00802)	0.008
Case 3	0.1	( $-\mathrm{0.0186}$, $\mathrm{0.0126}$)	(0.0077, 0.0082)
Case 4	0.5	( $-\mathrm{0.0914}$, $\mathrm{0.0632}$)	(0.0067, 0.0092)
Case 5	1	( $-\mathrm{0.1827}$, $\mathrm{0.1264}$)	(0.0056, 0.0105)
Case A	0	( $-\mathrm{3.2}$, $\mathrm{3.12}$ )	0.008	0.004
Case B	0	( $-\mathrm{1.56}$, $\mathrm{1.6}$)	0.008	0.010

Figure 8

Channel bottom offset ($h_{\mathrm{R}}$) and physical thalweg slope ($S_{\mathrm{0}}$) for synthetic test cases 2–5 that are random perturbations of the baseline smooth $S_{\mathrm{0}}=\mathrm{0.008}$ of Case 1.

[Figure omitted. See PDF]

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 20$$z_{\mathrm{0}}^{\left[c\right]}\left(x\right)=z_{\mathrm{0}}^{\left[\mathrm{1}\right]}\left(x\right)+{\mathit{\alpha }}^{\left[c\right]}H\left(x\right):c\in \left\{\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}\right\},$$ where $H\left(x\right)$ is a set of random-generated numbers within the range of $-\mathrm{0.126}\le H\left(x\right)\le \mathrm{0.183}$. The upper and lower limits of $H\left(x\right)$ were selected to prevent the occurrence of a locally adverse slope – such conditions can be handled by SPRNT-RS but can cause convergence problems for some models. The ${\mathit{\alpha }}^{\left[c\right]}$ 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., 21$$z_{\mathrm{R}}^{\left[c\right]}\left(x\right)=z_{\mathrm{0}}^{\left[\mathrm{1}\right]}\left(x\right):c\in \left\{\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}\right\}.$$ Thus, the SPRNT-RS simulations for cases 2–5 use uniform $S_{\mathrm{R}}$ over the reach such that the bed offset ($h_{\mathrm{R}}$) represents the physical geometric perturbations. Noting from Eq. () that $h_{\mathrm{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 22$$h_{\mathrm{R}}^{\left[c\right]}\left(x\right)={\mathit{\alpha }}^{\left[c\right]}H\left(x\right).$$

For synthetic test cases A and B in Table , 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_{\mathrm{0}}$) when $S_{\mathrm{R}}$ is set to an arbitrary value.

Forcing for all seven test cases is a constant inflow boundary of 283 m${}^{\mathrm{3}}$ s${}^{-\mathrm{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_{\mathrm{0}}=\mathrm{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${}^{\mathrm{2}}$ with total length of 10.7 km for the area illustrated in Fig. . Bathymetric survey data are available courtesy of the city of Austin (Fig. ). 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 cross-section spacing. Arguably, these issues are related to sharp changes in $A$ that lead to non-smooth source terms despite the RS method; however, this issue requires further investigation. Similar numerical instability behavior can also be found in the HEC-RAS 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 HEC-RAS 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 cross-section 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 step-pool transcritical elements in the downstream channel in which the HEC-RAS solution is strongly influenced by the ad hoc local partial inertia (LPI) algorithm . The smoothing introduced by LPI makes it difficult to draw conclusions from a comparison between SPRNT-RS and HEC-RAS across transcritical locations.

Figure 9

Main stem of Waller Creek and catchment in Austin (Texas, USA). Basemap © by Open Geospatial Consortium (OGC) Web Mapping Service (WMS).

[Figure omitted. See PDF]

Figure 10

Surveyed cross sections of Waller Creek. Only 140 out of 327 cross sections are shown for clarity. Elevations are relative to mean sea level. Data courtesy of the city of Austin.

[Figure omitted. See PDF]

A time-invariant upstream inflow boundary condition is set to 25 m${}^{\mathrm{3}}$ s${}^{-\mathrm{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_{\mathrm{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${}^{\mathrm{3}}$ s${}^{-\mathrm{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 SPRNT-RS and HEC-RAS models.

The thalweg elevation, $z_{\mathrm{0}}\left(x\right)$, and the reference elevation, $z_{\mathrm{R}}\left(x\right)$, of the RS method (as determined by the approximate spline fit described above) for Waller Creek are shown in Fig. a. The $z_{\mathrm{0}}$ and $z_{\mathrm{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 B-spline technique for generating $z_{\mathrm{R}}\left(x\right)$ does not force the overall reach slope to be identical. Because $z_{\mathrm{R}}\left(x\right)$ is mathematically arbitrary, there is no need to force an exact match. Although $z_{\mathrm{0}}\left(x\right)$ and $z_{\mathrm{R}}\left(x\right)$ are similar in Fig. a, $S_{\mathrm{0}}\left(x\right)$ from the raw data is discontinuous and varies over a wide range (up to 4 times the reach slope), as illustrated in Fig. b. Note that $S_{\mathrm{0}}\left(x\right)$ also includes negative slopes (i.e., adverse gradient sections) which can cause convergence problems for some numerical solvers. In contrast, as shown in Fig. c, the approximate cubic B-spline used to generate $S_{\mathrm{R}}\left(x\right)$ from $z_{\mathrm{R}}\left(x\right)$ 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 . The Lipschitz smoothness of $S_{\mathrm{0}}$ versus $S_{\mathrm{R}}$ can be better understood by evaluating the gradient of the slope, i.e., the second derivative of $z_{\mathrm{0}}$ and $z_{\mathrm{R}}$, as shown in Fig. d. The $S_{\mathrm{0}}$ formulation clearly lacks smoothness in the higher derivative.

Figure 11

Waller Creek bottom elevation and slope: (a) $z$ elevation, with 5 m subtracted from $z_{\mathrm{R}}$ for clarity; (b) $S_{\mathrm{0}}$ thalweg slope between cross sections; (c) $S_{\mathrm{R}}$ smoothed bottom slope; and (d) second derivative of $z_{\mathrm{0}}$ and $z_{\mathrm{R}}$. Note the $y$-axis scaling in (c) has reduced limits compared to (b) to better show the smoothness achieved by the spline fit.

[Figure omitted. See PDF]

Table 3

Data for model setup of Waller Creek test cases.

Case	Slope formulation	Slope range	Model usage
WC${}_{\mathrm{RS}}$	RS method	$\mathrm{0.0033}<S_{\mathrm{R}}<\mathrm{0.0147}$	SPRNT
WC${}_{\text{HEC-S}}$	Conventional	$-\mathrm{0.0328}<S_{\mathrm{0}}<\mathrm{0.0393}$	HEC-RAS (steady)
$W{C}_{\text{HEC-U}}$	Conventional	$-\mathrm{0.0328}<S_{\mathrm{0}}<\mathrm{0.0393}$	HEC-RAS (unsteady)

3.8 Analysis methods

To evaluate the performance of the RS method relative to conventional formulations, four depth-based indicators are employed, as described below. For these definitions, the control (superscript $\left[C\right]$) is the solution for the analytical test case and HEC-RAS results for the synthetic channel and Waller Creek test cases. Note that the synthetic tests use unsteady HEC-RAS, whereas the Waller Creek study uses comparisons to both steady and unsteady versions of the model. The test case (superscript $\left[T\right]$) is always the SPRNT-RS simulation. These measures can be considered error metrics for the comparison to the analytical solutions and difference metrics for model–model comparisons.

• Normalized difference ($\mathit{\rho }$). A nondimensional index to describe the local difference in depth can be defined as

  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)}},$$ 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 steady-state 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 ($\mathit{\zeta }$). The mean of the absolute value of $\mathit{\rho }\left(x\right)$ over the domain provides an overall nondimensional indicator of the depth error: 24$${\mathit{\zeta }}^{[T:C]}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{x=\mathrm{1}}^N\mid {\mathit{\rho }}^{[T:C]}\left(x\right)\mid ,$$ where $N$ is the total number of cross sections. We use the absolute value so that positive errors do not cancel negative errors, and $\mathit{\zeta }$ is a representative scale of the discrepancy between models.

• Mean absolute error (MAE). The overall dimensional error is characterized as 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 ,$$ and the nondimensional form of overall error is 26$$\begin{array}{rl}\text{MAE}& \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}$$

• Root mean square error (RMSE). A standard dimensional measure of the squared error is 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}}}.$$ The nondimensional form of RMSE is computed by the following equation: 28$$\begin{array}{rl}\text{RMSE}& \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}$$ Both the MAE and RMSE are also reported in nondimensional form where the normalization scale is presented.

4 Results

4.1 Analytical test cases

The water surface elevations for the analytical solutions and SPRNT-RS simulations are shown in Fig. . Visually, the analytical and simulated results across all six cases are identical. Error metrics following Sect.  are provided in Table . The normalized differences $\left(\mathit{\rho }\right)$ are less than 1 % and are consistent with absolute errors of $O\left({\mathrm{10}}^{-\mathrm{3}}\right)$ m, which are negligible compared with the water depths $\ge \mathrm{1}$ m. The spatial distributions of the normalized error are shown in Fig. . By comparing this figure with Fig. , it can be seen that $\mathit{\rho }\left(x\right)$ fluctuates with the change in bed slope. Similar behavior can also be found for model results reported in .

Figure 12

Comparison between simulated water surface elevation from SPRNT-RS and analytical solution for analytical test cases of .

[Figure omitted. See PDF]

Figure 13

Spatial distribution of normalized difference ($\mathit{\rho }$) for analytical test cases of .

[Figure omitted. See PDF]

Table 4

Difference measures using Eqs. ()–() for analytical test cases of . Nondimensional MAE and RMSE are shown in parentheses.

Case	Min ($\mathit{\rho }$)	Max ($\mathit{\rho }$)	$\mathit{\zeta }$	MAE (m)	RMSE (m)
UR1	$-\mathrm{0.47}\mathit{\%}$	$\mathrm{0.33}\mathit{\%}$	$\mathrm{0.09}\mathit{\%}$	$\mathrm{0.00104}$ $\left(\mathrm{0.032}\mathit{\%}\right)$	$\mathrm{0.00131}$ $\left(\mathrm{0.048}\mathit{\%}\right)$
UR3	$-\mathrm{0.73}\mathit{\%}$	$\mathrm{0.55}\mathit{\%}$	$\mathrm{0.28}\mathit{\%}$	$\mathrm{0.00385}$ $\left(\mathrm{0.044}\mathit{\%}\right)$	$\mathrm{0.00457}$ $\left(\mathrm{0.058}\mathit{\%}\right)$
UT1	$-\mathrm{0.19}\mathit{\%}$	$\mathrm{0.33}\mathit{\%}$	$\mathrm{0.07}\mathit{\%}$	$\mathrm{0.00085}$ $\left(\mathrm{0.027}\mathit{\%}\right)$	$\mathrm{0.00108}$ $\left(\mathrm{0.042}\mathit{\%}\right)$
UT2	$-\mathrm{0.73}\mathit{\%}$	$\mathrm{0.62}\mathit{\%}$	$\mathrm{0.28}\mathit{\%}$	$\mathrm{0.00388}$ $\left(\mathrm{0.064}\mathit{\%}\right)$	$\mathrm{0.00458}$ $\left(\mathrm{0.084}\mathit{\%}\right)$
VR1	$-\mathrm{0.13}\mathit{\%}$	$\mathrm{0.08}\mathit{\%}$	$\mathrm{0.03}\mathit{\%}$	$\mathrm{0.00055}$ $\left(\mathrm{0.006}\mathit{\%}\right)$	$\mathrm{0.00070}$ $\left(\mathrm{0.009}\mathit{\%}\right)$
VR2	$-\mathrm{0.15}\mathit{\%}$	$\mathrm{0.13}\mathit{\%}$	$\mathrm{0.06}\mathit{\%}$	$\mathrm{0.00090}$ $\left(\mathrm{0.026}\mathit{\%}\right)$	$\mathrm{0.00108}$ $\left(\mathrm{0.031}\mathit{\%}\right)$

4.2 Synthetic test cases

Results for the baseline synthetic test, Case 1, are shown in Fig. . The SPRNT-RS method produces visually the same solution to HEC-RAS with $z_{\mathrm{R}}=z_{\mathrm{0}}$. Similarly, the comparison of model results for depth ($h_{\mathrm{0}}$) for test cases 2–5 are visually indistinguishable, as shown in the left column of Fig. . The quantitative difference measures for the synthetic tests are provided in Table , and the spatial distributions of $\mathit{\rho }\left(x\right)$ are illustrated in the right column of Fig. . Values for $\mathit{\rho }\left(x\right)$ in cases 2 and 3 are slightly below zero ($\approx \mathrm{0.02}\mathit{\%}$) over the entire domain, indicating the SPRNT-RS solution has a slightly higher water surface than the HEC-RAS solution for small perturbations in the bed slope. With the increased bottom perturbations in cases 4 and 5, the $\mathit{\rho }\left(x\right)$ range is larger (and includes a positive range), but the bounding values are still trivial. The $\mathit{\zeta }$ and RMSE measures show that the nondimensional and dimensional overall differences are small. The MAE and RMSE climb slightly with the increasing $h_{\mathrm{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 ), and further results are not shown as they are visually identical to those for baseline Case 1 illustrated in Fig. . Note that the RMSE for both cases is identical to Case 1, which indicates the solution for the SPRNT-RS method with $S_{\mathrm{R}}\ne S_{\mathrm{0}}$ is very close to the baseline solution with $S_{\mathrm{R}}=S_{\mathrm{0}}$.

Figure 14

Simulated profile of water surface elevation (upper line) and channel bottom (lower line) for synthetic test Case 1 using SPRNT-RS (red) and unsteady HEC-RAS (black).

[Figure omitted. See PDF]

Figure 15

Water depth, $h_{\mathrm{0}}$, (left column) and normalized difference, $\mathit{\rho }$, (right column) for synthetic test cases with perturbed bathymetry.

[Figure omitted. See PDF]

Table 5

Difference measures between SPRNT-RS and HEC-RAS using Eqs. ()–() for synthetic test cases. Nondimensional MAE and RMSE are shown in parentheses.

Case	Min ($\mathit{\rho }$)	Max ($\mathit{\rho }$)	$\mathit{\zeta }$	MAE (m)	RMSE (m)
Case 1 (baseline)	$-\mathrm{0.107}\mathit{\%}$	$\mathrm{0.005}\mathit{\%}$	$\mathrm{0.0281}\mathit{\%}$	$\mathrm{0.00139}$ $\left(\mathrm{0.028}\mathit{\%}\right)$	$\mathrm{0.00176}$ $\left(\mathrm{0.035}\mathit{\%}\right)$
Case 2	$-\mathrm{0.031}\mathit{\%}$	$-\mathrm{0.012}\mathit{\%}$	$\mathrm{0.0251}\mathit{\%}$	$\mathrm{0.00124}$ $\left(\mathrm{0.025}\mathit{\%}\right)$	$\mathrm{0.00126}$ $\left(\mathrm{0.025}\mathit{\%}\right)$
Case 3	$-\mathrm{0.037}\mathit{\%}$	$-\mathrm{0.005}\mathit{\%}$	$\mathrm{0.0252}\mathit{\%}$	$\mathrm{0.00125}$ $\left(\mathrm{0.025}\mathit{\%}\right)$	$\mathrm{0.00129}$ $\left(\mathrm{0.026}\mathit{\%}\right)$
Case 4	$-\mathrm{0.073}\mathit{\%}$	$\mathrm{0.032}\mathit{\%}$	$\mathrm{0.0284}\mathit{\%}$	$\mathrm{0.00141}$ $\left(\mathrm{0.028}\mathit{\%}\right)$	$\mathrm{0.00168}$ $\left(\mathrm{0.034}\mathit{\%}\right)$
Case 5	$-\mathrm{0.107}\mathit{\%}$	$\mathrm{0.095}\mathit{\%}$	$\mathrm{0.0426}\mathit{\%}$	$\mathrm{0.00212}$ $\left(\mathrm{0.043}\mathit{\%}\right)$	$\mathrm{0.00249}$ $\left(\mathrm{0.049}\mathit{\%}\right)$
Case A	$-\mathrm{0.107}\mathit{\%}$	$\mathrm{0.005}\mathit{\%}$	$\mathrm{0.0281}\mathit{\%}$	$\mathrm{0.00139}$ $\left(\mathrm{0.028}\mathit{\%}\right)$	$\mathrm{0.00175}$ $\left(\mathrm{0.035}\mathit{\%}\right)$
Case B	$-\mathrm{0.107}\mathit{\%}$	$\mathrm{0.005}\mathit{\%}$	$\mathrm{0.0281}\mathit{\%}$	$\mathrm{0.00139}$ $\left(\mathrm{0.028}\mathit{\%}\right)$	$\mathrm{0.00175}$ $\left(\mathrm{0.035}\mathit{\%}\right)$

4.3 Waller Creek test case

Waller Creek has been simulated with SPRNT-RS (denoted as WC${}_{\text{RS}}$ in the following figures), the HEC-RAS unsteady solver (WC${}_{\text{HEC-U}}$), and the HEC-RAS steady solver (WC${}_{\text{HEC-S}}$). Figure  shows water surface elevations for SPRNT-RS and unsteady HEC-RAS. For clarity, the upper 40% of the domain (which has similar good behavior) is not shown. Figure  shows the spatial distribution of the normalized difference, $\mathit{\rho }$, for these simulations. The differences are roughly within $\pm \mathrm{4}\mathit{\%}$ across the entire domain. The maximum and minimum difference both occur at two adjacent nodes close to 7800 m with $\mathrm{4.14}\mathit{\%}$ and $-\mathrm{3.07}\mathit{\%}$, respectively. Figure  provides a similar comparison of water surface elevations between SPRNT-RS and the steady HEC-RAS case. The results are visually quite similar to the comparison with unsteady HEC-RAS. A direct comparison of surface elevations for unsteady and steady HEC-RAS does not provide any further insight and is omitted for brevity. However, to quantitatively evaluate the differences between SPRNT-RS and HEC-RAS, it is useful to compute difference measures between the unsteady and steady HEC-RAS models themselves, as well as the differences between SPRNT-RS and both models, as provided in Table . Overall, the SPRNT-RS result have marginally better consistency with unsteady HEC-RAS than with steady HEC-RAS. Of greater importance is that the behavior of SPRNT-RS relative to unsteady HEC-RAS has the same order of differences as the comparison of unsteady HEC-RAS to steady HEC-RAS. These results imply that the differences between SPRNT-RS and unsteady HEC-RAS are reasonable for the different numerical methods given the geometric variability of Waller Creek.

Figure 16

Comparison of SPRNT-RS to unsteady HEC-RAS for water surface elevations in Waller Creek simulations with expanded detail to show similarities. For clarity, 5 m is subtracted from the channel $z$ elevations.

[Figure omitted. See PDF]

Figure 17

Normalized difference $\mathit{\rho }\left(x\right)$ between SPRNT-RS and unsteady HEC-RAS.

[Figure omitted. See PDF]

Figure 18

Comparison of SPRNT-RS to steady HEC-RAS for water surface elevations in Waller Creek simulations with expanded detail to show similarities. For clarity, 5 m is subtracted from the channel $z$ elevations.

[Figure omitted. See PDF]

Table 6

Difference metrics for Waller Creek simulation results. Nondimensional MAE and RMSE are shown in parentheses.

Case comparison	Min ($\mathit{\rho }$)	Max ($\mathit{\rho }$)	$\mathit{\zeta }$	MAE (m)	RMSE (m)
WC${}_{\text{RS}}$ : WC${}_{\text{HEC-U}}$	$-\mathrm{3.07}$ %	4.14 %	$\mathrm{0.85}\mathit{\%}$	$\mathrm{0.056}$ $\left(\mathrm{0.883}\mathit{\%}\right)$	$\mathrm{0.077}$ $\left(\mathrm{1.284}\mathit{\%}\right)$
WC${}_{\text{RS}}$ : WC${}_{\text{HEC-S}}$	$-\mathrm{8.43}$ %	4.58 %	$\mathrm{1.25}\mathit{\%}$	$\mathrm{0.081}$ $\left(\mathrm{1.287}\mathit{\%}\right)$	$\mathrm{0.122}$ $\left(\mathrm{1.990}\mathit{\%}\right)$
WC${}_{\text{HEC-U}}$ : WC${}_{\text{HEC-S}}$	$-\mathrm{2.70}$ %	7.03 %	$\mathrm{1.30}\mathit{\%}$	$\mathrm{0.086}$ $\left(\mathrm{1.301}\mathit{\%}\right)$	$\mathrm{0.128}$ $\left(\mathrm{1.931}\mathit{\%}\right)$

5 Discussion

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_{\mathrm{0}}$) as a control model (see Sect. ). Invariably, discontinuous topography for SPRNT without RS caused either an oscillatory solution or numerical instability. In contrast, both HEC-RAS (using $\mathit{\eta }$) and SPRNT-RS (using $S_{\mathrm{R}}$) provide stable, non-oscillatory solutions.

The analytical results in Sect. , supplemented by additional results in , validate the SPRNT-RS 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_{\mathrm{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 . These results imply that the transformation from the conventional $h_{\mathrm{0}}$ and $S_{\mathrm{0}}$ form of the SVEs to the $h_{\mathrm{a}}$ and $S_{\mathrm{R}}$ form of Eq. () 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.  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_{\mathrm{R}}$. The arbitrary selection of an $S_{\mathrm{R}}\ne S_{\mathrm{0}}$ results in identical solutions to $S_{\mathrm{R}}=S_{\mathrm{0}}$. Secondly, the synthetic test cases show that the SPRNT-RS method can be applied with non-smooth geometry at the discretization scale (i.e., random perturbations of the physical bottom slope), which caused non-convergent behavior in the baseline SPRNT model. As a control, we have compared SPRNT-RS with the accepted HEC-RAS model that remains stable for these test cases as it is solved with the piezometric pressure gradient rather than splitting into $S_{\mathrm{0}}$ and the gradient of $h_{\mathrm{0}}$. The results indicate that SPRNT-RS provides numerical solutions that are nearly identical to HEC-RAS for the non-smooth geometry test cases. Thus, using a Lipschitz smooth $S_{\mathrm{R}}$ provides a stable numerical solution for non-smooth geometry without altering the physical representation of non-smooth geometry.

The Waller Creek test case in Sect.  provides a more challenging comparison of SPRNT-RS to HEC-RAS. For this test case, the geometry discontinuities include adverse slopes and local $S_{\mathrm{0}}$ that are $\pm \mathrm{400}\mathit{\%}$ of the reach-average slope, which contrasts with perturbations of $\pm \mathrm{30}\mathit{\%}$ used in the synthetic test cases of Sect. . Again, SPRNT-RS is shown to be close to the unsteady HEC-RAS solution. The model differences are within reasonable ranges, as illustrated by the fact that they are similar to the differences between HEC-RAS steady and unsteady versions. Nevertheless, it remains possible that the minor differences between HEC-RAS and SPRNT-RS 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 HEC-RAS that are not present in SPRNT-RS. Specifically, notes that for computational efficiency and to reduce “troublesome convergence problems at discontinuities in the river geometry,” the unsteady HEC-RAS code uses a linearization technique developed by and – note that the latter document is cited by 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 HEC-RAS. Unfortunately, we do not have direct access to the HEC-RAS code and thus rely on the discussion of HEC-RAS stability in the literature and the methodology in HEC-RAS manuals .

5.2

Why not just use $\partial \mathit{\eta }/\partial x$?

One might wonder whether $S_{\mathrm{R}}$ or $S_{\mathrm{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_{\mathrm{R}}$, it is worth considering why $S_{\mathrm{0}}$ is presently used. We have not been able to determine exactly when $S_{\mathrm{0}}$ was first used with the SVEs, but from a hydrology viewpoint $S_{\mathrm{0}}$ provides consistency between kinematic wave solutions (which use $S_{\mathrm{0}}=S_{\mathrm{f}}$) and the SVEs. Thus including $S_{\mathrm{0}}$ is a logical step when considering reduced-physics approaches . Arguably, a well-chosen $S_{\mathrm{R}}$ that matches the large-scale $S_{\mathrm{0}}$ will serve the same purpose. The $S_{\mathrm{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.,. For these model, the advantage of the $S_{\mathrm{0}}$ form is that, when $S_{\mathrm{0}}$ equals 0 and $S_{\mathrm{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_{\mathrm{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 $\mathit{\eta }$ approach is obvious, it has a key disadvantage when applied in large-scale simulations. Over large distances the free surface is monotonically increasing upstream, which has consequences for employing implicit or semi-implicit numerical solutions in a continental river dynamics framework . Briefly, when modeling a river system from an estuary ($\mathit{\eta }\sim \mathrm{0}$ m) to mountain headwaters $(\mathit{\eta }\sim {\mathrm{10}}^{\mathrm{3}}$ m), the solution variable, $\mathit{\eta }$, nominally covers 3 orders of magnitude. Furthermore, as local variations on the order of ${\mathrm{10}}^{-\mathrm{2}}$ m affect the hydrostatic pressure gradient, the solution of $\mathit{\eta }$ 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 down-slope body force ($S_{\mathrm{R}}$) and a local residual ($\partial h_{\mathrm{a}}/\partial x)$ is effectively removing a large-scale gradient from the solution variables, which will generally improve numerical behavior.

Despite the above disadvantages, the $\mathit{\eta }$ form retains some advantages in creating conservative finite-volume formulations of the Saint-Venant equations . Arguably, such methods should be confined to explicit time-marching schemes or localized solutions, in which $\mathit{\eta }$ 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 SPRNT-RS approach and most, if not all, conventional models (including unsteady HEC-RAS) is that our method algebraically revises the Saint-Venant 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/SPRNT-RS) and only a single model (HEC-RAS) 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 well-balanced 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 well-balanced 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 , which can be suppressed with the ad hoc local partial inertia (LPI) scheme of . 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_{\mathrm{0}}$ to $S_{\mathrm{R}}$ is algebraically exact, the application of the RS method requires some method to select the distributed $z_{\mathrm{R}}\left(x\right)$ and to determine $S_{\mathrm{R}}\left(x\right)$. Poor selection of $z_{\mathrm{R}}$ can theoretically result in non-smooth $S_{\mathrm{R}}$. In the present work, the profile of $z_{\mathrm{R}}$ in the Waller Creek case is generated by the cubic B-spline 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_{\mathrm{R}}\left(x\right)$. Our results indicate that approximating cubic B-splines is adequate for producing smooth $z_{\mathrm{R}}$ for the tested geometries and that the solutions are robust for the selection of $z_{\mathrm{R}}$ as long as $S_{\mathrm{R}}$ is smooth . However, it is likely there are limitations to applying the RS method in large-scale river network simulations that will make it difficult to use a simple globally applied knot spacing. Such networks might consist of ${\mathrm{10}}^{\mathrm{4}}$ to ${\mathrm{10}}^{\mathrm{5}}$ reaches spanning wide geographical regions with a variety of topology and inconsistent data availability. Some reaches may have well-defined cross sections at close spacing, other reaches might be poorly documented . 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_{\mathrm{R}}\left(x\right)$, 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. . The baseline governing equations for SPRNT are of the form of Eq. (), the so-called “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 , 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. () can arguably be applied in the Cunge–Liggett form or any other conservative form of the SVEs. Similarly, the fundamental algebraic transformation to $S_{\mathrm{R}}$ and $\partial h_{\mathrm{a}}/\partial x$ will be equally valid in any finite-volume method using $S_{\mathrm{0}}$ and $h_{\mathrm{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_{\mathrm{a}}$, $h_{\mathrm{R}}$, and $z_{\mathrm{R}}$ in place of $h_{\mathrm{0}}$ and $z_{\mathrm{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 well-balanced problem for conservative forms as reviewed in and . 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.

6 Conclusions

The reference slope (RS) method introduces a new form of the Saint-Venant 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_{\mathrm{R}}\left(x\right)$, with computed reference slopes, $S_{\mathrm{R}}\left(x\right)$, and associated depths, $h_{\mathrm{a}}\left(x\right)$. 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 open-source Saint-Venant solver as SPRNT-RS. In this study, SPRNT-RS was compared to both analytical solutions and the conventional HEC-RAS 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 HEC-RAS that do not appear in SPRNT-RS.

As discussed in the Sect. , when faced with non-smooth 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 cross-section 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 fine-resolution reach-scale modeling to large-scale continental river dynamic simulations or develop massively parallel stormwater network models for megacities

The RS method is not specific to SPRNT but can be adapted to any Saint-Venant solver that uses a bottom slope ($S_{\mathrm{0}}$) term in the discretization. The mathematical change is conceptually trivial, but the actual effort depends on how cross-section geometry is embedded in the code. The code for both SPRNT and SPRNT-RS are available under open-source license at GitHub .

Appendix A Geometry adjustments for stable unsteady HEC-RAS simulations

As discussed in Sect. , the stability of the SPRNT-RS simulations for Waller Creek was ensured by merging 36 computational elements in which the cross-section spacing was closer than 10 m. This minimum spacing cutoff was selected as being well below the median spacing of $\mathrm{28.6}$ m and mean spacing of $\mathrm{37.7}$ m and proved adequate for ensuring SPRNT-RS stability in the tested simulations. Unfortunately, stability of unsteady HEC-RAS required the further removal of three cross-section elements (shown in Fig. ) and reducing Manning's $n$ at six additional cross sections (listed in Table ). Selecting these changes was a matter of art rather than science as we could not identify clear criteria for cross-section removal or Manning's $n$ adjustments for HEC-RAS – other than that these locations appeared to be where instabilities appeared in unsteady HEC-RAS model runs. Although SPRNT-RS could run without these changes, for consistency in the model comparisons the geometry of the SPRNT-RS model was modified to exactly match the adjusted geometry required for the HEC-RAS model.

Figure A1

Cross sections removed from Waller Creek data set to provide numerical stability of unsteady HEC-RAS.

[Figure omitted. See PDF]

Table A1

Modified Manning's $n$ for cross-section stations in Waller Creek data set to provide numerical stability of unsteady HEC-RAS.

Station	Reach location	Original	Modified
number	(m)	$n$	$n$
30104	1384	0.06	0.04
30014	1412	0.06	0.04
29871	1454	0.055	0.04
29752	1490	0.06	0.04
29647	1522	0.055	0.04
29482	1572	0.055	0.04

Appendix B Notation

$A$	cross-section area (m${}^{\mathrm{2}}$)
$g$	gravitational acceleration (ms${}^{-\mathrm{2}}$)
$h_{\mathrm{0}}$	water depth (m)
$h_{\mathrm{a}}$	associated water depth (m)
$h_{\mathrm{R}}$	reference height (m)
$n$	Gauckler–Manning–Strickler roughness (${\text{m}}^{-\mathrm{1}/\mathrm{3}}\text{s}$)
$P_{\mathrm{w}}$	wetted perimeter (m)
$Q$	volumetric flow rate (${\text{m}}^{\mathrm{3}}{\text{s}}^{-\mathrm{1}}$)
$q_{\mathrm{\ell }}$	flow rate per unit length through channel sides (m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$)
$S_{\mathrm{0}}$	channel bottom slope
$S_{\mathrm{f}}$	channel friction slope
$S_{\mathrm{R}}$	channel reference slope
$S_{\mathrm{SW}}$	channel sidewall slope
$t$	time (s)
$W$	channel width (m)
$W_{\mathrm{B}}$	channel bottom width ($\text{(}m)$)
$x$	along-channel spatial coordinate
$z_{\mathrm{0}}$	channel bottom elevation (m)
$z_{\mathrm{R}}$	reference elevation (m)
$\mathit{\alpha }$	bottom displacement coefficient
$\mathit{ϵ}$	second derivative of water surface elevation (m${}^{-\mathrm{1}}$)
$\mathit{\eta }$	water surface elevation (m)
$\mathit{\rho }$	normalized difference between results
$\mathit{\zeta }$	absolute mean normalized difference (AMND)

Code availability

The complete code for the reference slope module and SPRNT is available at Github (https://github.com/frank-y-liu/SPRNT, 10.5281/zenodo.3978737, ).

Data availability

All test case files and results have been uploaded to a public repository under the Texas Data Repository and can be accessed through 10.18738/T8/BXJBF5.

Author contributions

CWY and BRH designed and performed the experiments, analyzed the data, and wrote the paper. FL wrote the code for the reference slope method.

Competing interests

The authors declare that they have no conflict of interest.

Disclaimer

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.

Acknowledgements

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 UT-Battelle, LLC, for the U.S. Department of Energy.

Financial support

This research has been supported by the U.S. National Science Foundation (grant no. CCF-1331610) and the U.S. Environmental Protection Agency (cooperative agreement no. 83595001).

Review statement

This paper was edited by Roberto Greco and reviewed by two anonymous referees.