Conservative finite-volume forms of the Saint-Venant equations for hydrology and urban drainage

New integral, finite-volume forms of the Saint-Venant equations for one-dimensional (1-D) open-channel flow are derived. The new equations are in the flux-gradient conservation form and transfer portions of both the hydrostatic pressure force and the gravitational force from the source term to the conservative flux term. This approach prevents irregular channel topography from creating an inherently non-smooth source term for momentum. The derivation introduces an analytical approximation of the free surface across a finite-volume element (e.g., linear, parabolic) with a weighting function for quadrature with bottom topography. This new free-surface/topography approach provides a single term that approximates the integrated piezometric pressure over a control volume that can be split between the source and the conservative flux terms without introducing new variables within the discretization. The resulting conservative finite-volume equations are written entirely in terms of flow rates, cross-sectional areas, and water surface elevations – without using the bottom slope (S₀). The new Saint-Venant equation form is (1) inherently conservative, as compared to non-conservative finite-difference forms, and (2) inherently well-balanced for irregular topography, as compared to conservative finite-volume forms using the Cunge–Liggett approach that rely on two integrations of topography. It is likely that this new equation form will be more tractable for large-scale simulations of river networks and urban drainage systems with highly variable topography as it ensures the inhomogeneous source term of the momentum conservation equation is Lipschitz smooth as long as the solution variables are smooth.

1 Introduction

The one-dimensional (1-D) Saint-Venant equations (SVEs) are the simplest equations that capture the full dynamics of river and open-channel flow, and yet they are not universally used where river networks are explicitly represented in hydrological and urban drainage models. Furthermore, where the SVEs are used the numerical solution methods typically employ a non-conservative form, which inherently has greater numerical dissipation than the conservative form. But even the undesirable non-conservative form of the SVEs is rarely used for large systems – instead, reduced-physics models are common. These approaches a priori neglect some part of the flow dynamics for simpler computation (see Sects. 2 and 3 for discussions of these issues). In their call for “hyperresolution global land surface modeling”, Wood et al. (2011) postulated that the principal issues that limit SVE use in river networks could be resolved by application of more computing power and more precise river geometry. Hodges (2013) argued that there were wider challenges to using the SVEs, but these are not sufficient reasons relying on reduced-physics models that are calibrated to get the right answer for the wrong reason. In the intervening years some progress has been made toward improving applications of the SVEs, but large-scale hydrological modeling of river networks with application of supercomputer power continues to rely on reduced-physics methods (e.g., the US National Water Model, Cohen et al., 2018). Arguably, the dynamics of river flow are the most observable and should be the easiest part of hydrology to model (or at least with lower uncertainty in boundaries and forcing compared to overland or groundwater flow), so it is unsatisfactory that more than a half century after Preissmann (1961) we still are not universally using the SVEs in hydrological modeling of river networks.

The present work evolved out of a frustration with the slow pace of improvement in SVE modeling. Taking a step backwards, we can ask the following: is there something fundamental in the common forms of the SVEs that hinders progress? Motivated by an analysis of the equation forms (Sect. 2) and a study of the wealth of past work in the SVEs (Sect. 3), new insights were developed and are presented herein. The fundamental theses of the present work are as follows: (1) conservative formulation of equations should be used for the next generation of river network models, and (2) the appearance of the channel slope (S₀) in the SVEs for channels with irregular topography is a principal cause of instabilities and extended computational time. Neither thesis can be demonstrated herein – this work is merely a first step that provides the theoretical foundations for a conservative and inherently well-balanced approach that highlights the minimal level of approximations needed for a SVE form with irregular topography. It remains for future studies to compare models built on these foundations to the existing approaches to determine if the new forms provide significant numerical advantages.

The new conservative form of the SVEs is developed with a goal of addressing challenges associated with modeling large-scale 1-D flow network systems. In the process of developing the new form, we will encounter a philosophical question as to whether the primary vertical variable in a large-scale network solution should be the depth (H) or the water surface elevation (η). Despite this author's prior work with H primacy (Liu and Hodges, 2014), we shall see that there are advantages to using η, which is identical to the piezometric pressure and hence uniform over a channel cross-sectional area. The quadrature of the subgrid piezometric pressure gradient and subgrid-scale topography can be handled in a single new term that is derived herein. Through this term interesting possibilities for analytically including hyperresolution bathymetric knowledge while retaining larger computational elements for large-scale modeling arise. This idea is not fully exploited within the present work, but the framework is developed for others to build upon.

In the remainder of this paper, Sect. 2 provides motivation and context in the differential forms of the SVEs. Section 3 provides a further overview of SVE modeling in the wide range of conservative and non-conservative forms. A new and complete derivation of the finite-volume form of the conservative 1-D momentum equation with minimal approximations is provided in Sect. 4. Approximate forms of the 1-D SVEs are presented in Sect. 5 and the final form of the equations and a discussion of their potential use are provided Sect. 6.

2 Motivation

In reach-scale hydraulic studies, the Saint-Venant equations are almost always solved in a conservative form (e.g., Karelsky et al., 2000; Lai et al., 2002; Papanicolaou et al., 2004; Sanders, 2001) but usually in a non-conservative form when used in river network hydrology and urban drainage networks (e.g., Liu and Hodges, 2014; Pramanik et al., 2010; Rossman, 2017; Saleh et al., 2013). Arguably, the reasons are in the difficulty in obtaining a well-balanced model for the conservative form and the inherent complexities/uncertainties of channel geometry across a large network (see discussion in Sect. 3). In general, conservative equation forms are valued as they ensure (with careful discretization) that transport modeling does not numerically create or destroy the transported variable. Indeed, the use of the conservative form for mass conservation is universal in models from hydraulics to hydrology – it is only for momentum that the non-conservative form remains common.

To set the context for this paper, consider the non-conservative form of the momentum equation that has been used in large river network solutions (Liu and Hodges, 2014)

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\partial Q}{\partial \mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial x}}\left(\mathit{\beta }{\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)+gA{\displaystyle \frac{\partial H}{\partial x}}=gA\left(S_{\mathrm{0}}-S_{\mathrm{f}}\right),\end{array}$$

where Q is the flow rate, A is the channel cross-sectional area, β is the momentum coefficient (associated with nonuniform velocities integrated over A), g is gravity, H is the water depth, S₀ is the channel slope, and S_f is the friction slope. The above equation has immediate physical appeal as each term represents a clearly understood piece of physics; i.e., the rate of change of the flow is affected by the gradient of nonlinear advection, the hydrostatic pressure gradient (driven by water depth), the gravitational force along the slope, and frictional resistance. The equation includes a conservative form of nonlinear advection (all the variables inside the gradient), but $gA\partial H/\partial x$ is non-conservative (that is, A, which is a function of H, is outside the gradient). The terms to the right-hand side (RHS) of the equal sign are considered source and sink terms that reflect creation and destruction of momentum. Thus, a conservative form of the above could be formally written as

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{\partial Q}{\partial \mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial x}}\left(\mathit{\beta }{\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)=-gA{\displaystyle \frac{\partial H}{\partial x}}+gA\left(S_{\mathrm{0}}-S_{\mathrm{f}}\right),\end{array}$$

where the trivial act of moving the pressure gradient term to the RHS is a recognition that it can cause non-conservation (source and sink) of the conserved momentum fluxes from the left-hand side (LHS). That is, any component on the RHS is capable of creating or destroying momentum, whereas the terms on the LHS (if properly discretized) cannot. The above equations can be contrasted with an equivalent momentum equation using the free-surface elevation:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle \frac{\partial Q}{\partial \mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial x}}\left(\mathit{\beta }{\displaystyle \frac{Q^{\mathrm{2}}}{A}}\right)=-gA{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}-gA{S}_{\mathrm{f}},\end{array}$$

which is obtained by substitution of $\mathit{\eta }=H+z_{\mathrm{b}}$ with z_b as the bottom elevation and $S_{\mathrm{0}}=-\partial z_{\mathrm{b}}/\partial x$. The equation is again conservative by virtue of including the gradient of the free-surface elevation as a source term (Ying et al., 2004; Wu and Wang, 2007; Ying and Wang, 2008). Comparison of Eqs. (2) and (3) shows how the introduction of S₀ affects the equation form. Of particular note is that we expect S_f=f(Q, A), so the latter equation will have a smooth source term as long as the solution variables themselves are smooth. In contrast, the smoothness of the source term in Eq. (2) inherently depends on the smoothness of the product AS₀ and compensation by the solution in $A\partial H/\partial x$ and AS_f. The key point is that smoothness in the source term is a mathematical necessity for the numerical solution of a partial differential equation to be well-posed (Iserles, 1996), but introduction of S₀ can place smoothness at the mercy of how well the numerical scheme responds to non-smooth forcing.

In general, there is an advantage to having as much of the physics as possible included on the flux-conservative side of the equation, which helps reduce difficulties associated with discretizations of the source terms (e.g., Pu et al., 2012; Vazquez-Cendon, 1999). The standard form of the conservative SVEs is arguably the form provided by Cunge et al. (1980), based on a derivation of Liggett (1975):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle \frac{\partial Q}{\partial \mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial x}}\left(\mathit{\beta }{\displaystyle \frac{Q^{\mathrm{2}}}{A}}+g{I}_{\mathrm{1}}\right)=g{I}_{\mathrm{2}}+gA\left(S_{\mathrm{0}}-S_{\mathrm{f}}\right),\end{array}$$

where I₁ and I₂ are integrated hydrostatic pressure forces across the channel

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{I}_{\mathrm{1}}=\underset{H}{\int }(H-z)B\mathrm{d}z\end{array}$$

and along the channel gradients

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{I}_{\mathrm{2}}=\underset{H}{\int }(H-z){\displaystyle \frac{\partial B}{\partial x}}\mathrm{d}z,\end{array}$$

with H(x) as the water depth, B(x, z) as the channel breadth as a function of elevation and along-stream location, and z as a coordinate direction measured from a common horizontal baseline in a direction opposite to gravitational acceleration. This form is also used by others with slightly different nomenclature but the same integral terms (e.g., Hernandez-Duenas and Beljadid, 2016; Sanders, 2001; Saavedra et al., 2003). It is convenient herein to call this the Cunge–Liggett form of the SVEs. The key point for both these terms is that they measure the interaction between the free surface and the channel shape; e.g., I₁ could also be written as

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{I}_{\mathrm{1}}=\underset{z_{\mathrm{b}}}{\overset{\mathit{\eta }}{\int }}(\mathit{\eta }-z)B\mathrm{d}z,\end{array}$$

where z_b is the channel bottom.

By comparing Eq. (4) with Eq. (2), we see that the novelty of Cunge–Liggett is in moving a portion of the $gA\partial H/\partial x$ from the RHS source term into the conservative flux on the LHS. Indeed, with this idea we see that Eq. (2) can be used with the product rule of differentiation to generate a slightly different conservative form:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle \frac{\partial Q}{\partial \mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial x}}\left(\mathit{\beta }{\displaystyle \frac{Q^{\mathrm{2}}}{A}}+gAH\right)=gH{\displaystyle \frac{\partial A}{\partial x}}+gA\left(S_{\mathrm{0}}-S_{\mathrm{f}}\right).\end{array}$$

Similar to the Cunge–Liggett form, the above uses a mathematical trick to place one part of the hydrostatic pressure force within the conservative flux gradient, while retaining the remainder as a source term. Thus, we see that the Cunge–Liggett form is not canonical, nor is it a form that necessarily better represents the physics. It is merely a form that is (sometimes) convenient for splitting the gradient of the total hydrostatic pressure force into conservative flux and non-conservative source terms.

The above brings up a question: if it is good to shift a portion of the hydrostatic pressure from the source to the conservative term as in Eqs. (4) and (8), then why not some or all of the gravitational potential associated with gAS₀? Underlying the Cunge–Liggett form and much (but not all) of the literature is the idea that gAS₀ is a source term that creates and destroys momentum. But this is also true of the hydrostatic pressure gradient and yet we commonly treat a portion within the convective flux. Thus, if Cunge–Liggett Eq. (4) is preferred over the baseline conservative form of Eq. (2) because a portion of the hydrostatic pressure is moved from the source term to the conservative flux term, then an equation that moves some (or all) of the gravitational potential from the source to the flux term should be equally valued. If this argument is accepted, then the preferred differential form for natural channels with the SVEs is none of those presented above, but perhaps an equation of the general form

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle \frac{\partial Q}{\partial \mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial x}}\left\{\mathit{\beta }{\displaystyle \frac{Q^{\mathrm{2}}}{A}}+f_{\mathrm{1}}(A,\mathit{\eta })\right\}=f_{\mathrm{2}}(A,\mathit{\eta })-gA{S}_{\mathrm{f}},\end{array}$$

where f₁ and f₂ are some functions (as yet unspecified) of the free-surface elevation rather than the depth. The free surface (η) can be interpreted as the uniform piezometric pressure over the cross-sectional area A, so the f₁ and f₂ functions are a generic way of splitting the piezometric pressure gradient between the source term and the conservative flux term, similar to how Cunge–Liggett handles the hydrostatic pressure gradient. Note that this proposed general approach removes the need for specifying S₀, which is important both in terms of source-term smoothness and producing well-balanced methods as discussed in detail in Sect. 3, below.

In this paper, complete derivations are presented to show that finite-volume formulations of the SVEs can be generated that are consistent with the general conservative differential form of Eq. (9). The new derivation is intended to bridge the gap between approaches used in high-resolution hydraulic models and those used in large-scale hydrology and urban drainage networks. The derivation provides a form of the SVEs that has mathematical rigor while preserving the simplicity of the non-conservative finite-difference discretizations that are common in hydrological and urban drainage literature. Herein, we focus only on the detailed presentation of the new equation form, reserving demonstration in a numerical model to future papers.

3 Background

3.1 Origination and use of the SVEs

The original equations of de Saint-Venant (de Saint-Venant, 1871) were written as

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}\left(AU\right)}{\mathrm{d}x}}=\mathrm{0},\text{(11)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{d}\mathit{\eta }}{\mathrm{d}x}}={\displaystyle \frac{\mathrm{1}}{g}}{\displaystyle \frac{\mathrm{d}U}{\mathrm{d}t}}+{\displaystyle \frac{U}{g}}{\displaystyle \frac{\mathrm{d}U}{\mathrm{d}x}}+{\displaystyle \frac{{\mathrm{\ell }}_{\mathrm{p}}}{A}}{\displaystyle \frac{F}{\mathit{\rho }g}},\end{array}$$

where U is the velocity, ℓ_p is the wetted perimeter, and F is the frictional force per unit bottom area along the channel, and other terms are as previously noted. We have taken the liberty of replacing Saint-Venant's notation of w, ζ, χ, s with the more modern nomenclature of A, η, ℓ_p, x, but otherwise have retained the original form. The momentum equation of de Saint-Venant, Eq. (11), is identical to Eq. (1) if we use Q=AU, $S_{\mathrm{f}}={\mathrm{\ell }}_{\mathrm{p}}F/\left(g\mathit{\rho }A\right)$, integrate over a cross section with the β coefficient, apply some calculus with the continuity equation, and define $\mathit{\eta }\equiv H+z_{\mathrm{b}}$ along with $S_{\mathrm{0}}=-\partial z_{\mathrm{b}}/\partial x$. From a practical perspective, the only thing that a hydrologist really needs to change from the original equation set is to replace the zero on the right-hand side of Eq. (10) with a source term representing the inflow/outflow per unit length from/into the catchment and groundwater. However, from a numerical modeling perspective, Eq. (11) is fundamentally non-conservative and suitable only for discretization in finite-difference forms. Although the full equation set is sometimes called the SVEs, for convenience in exposition we will use SVE as a shorthand for the momentum equation alone.

The SVEs are ubiquitous in the literature for a wide range of work and have a foundational role in flow-routing schemes in hydrological models and channel network models for urban drainage. However, there is an interesting gap between the equation forms used in large-scale systems and those used in shorter single-reach studies or modeling hydraulic features. For computational simplicity, large-scale network flow models often use a reduced set of equations, such as the Muskingum, kinematic wave, or the local inertia form (e.g., Wang et al., 2006; David et al., 2011, 2013; Getirana et al., 2017). Herein, we will follow the arguments of Hodges (2013) that we should be using the full SVEs; i.e., reduced-physics models should be seen as a stopgap measure as we wrestle with obtaining satisfactory SVE solution methods. As computational power has increased, our large-scale models have been moving towards the full SVEs but typically in a non-conservative form (e.g., Paiva et al., 2013; Liu and Hodges, 2014). For urban drainage networks, the US EPA Storm Water Management Model (SWMM) and variants built on this public domain model use a non-conservative finite-difference form of the SVEs. This model is widely applied (e.g., Gulbaz and Kazezyilmaz-Alhan, 2013; Hsu et al., 2000; Krebs et al., 2013); however deficiencies in conservation are a recognized problem (Rossman, 2017) and the SVE solver is the critical computational expense in the modeling system (Burger et al., 2014). Engineering river hydraulics problems are often solved using the US Army Corps of Engineers HEC-RAS software, which has free model executables with a proprietary (closed-source) code base. HEC-RAS uses a non-conservative finite-difference form of the SVEs based on methods pioneered in last quarter of the 20th century (Brunner, 2010). In contrast, more recent research models of the SVEs at short river-reach scales have typically used the equations presented as hyperbolic conservation laws that ensure both conservation and well-behaved solutions for subcritical, supercritical, and transcritical flows (e.g., Guinot, 2009; Ivanova et al., 2017; Papanicolaou et al., 2004; Sanders et al., 2003).

There is also a vast gulf between the spatial discretization of SVEs for large systems and smaller system studies in the hydraulics and applied mathematics literature (although the gap is getting narrower). For example in 2003 the SVEs were solved at 1 to 4 km spacing for 156 km of river (Saavedra et al., 2003). Seven years later we find 4 km cross-section spacing for 5×10³ km of river (Pramanik et al., 2010). By 2014 the state of the art was 100 m spacing for 15×10³ km of river (Liu and Hodges, 2014). In contrast, hydraulic studies have typically focused on 1 to 10 m spacing for 1 to 2 km test cases (e.g., Gottardi and Venutelli, 2003; Kesserwani et al., 2009; Sart et al., 2010; Venutelli, 2006). Between these extremes, single-reach river models with natural geometry are typically modeled over river lengths less than 20 km with grid cells on the order of 10 m to more than 100 m (Sanders et al., 2003; Catella et al., 2008; Castellarin et al., 2009; Lai and Khan, 2012).

3.2 Preissmann vs. Godunov

Computational modeling of the SVEs is arguably a long-running contest between the differential, finite-difference governing equations pioneered by Preissmann (1961) and the integral, finite-volume formulations derived from Godunov (1959). Clearly this is a simplification as there are wide-ranging contributions across both numerical methods and implementation schemes – but the literature is simply too broad to discuss all developments in anything less than a book. Nevertheless, using a dialectic of Preissmann vs. Godunov is a useful way of thinking about the major developments and providing context for the equations derived herein. Preissmann developed effective finite-difference methods applied to the non-conservative form of the equations, which is the first basis for comparison of succeeding finite-difference models. The introduction of Roe's approximate Riemann solver (Roe, 1981) and analyses of Harten et al. (1983) made Godunov-like methods tractable and set off a multi-decadal development of finite-volume methods. The literature with these two methods is vast, but a reasonable cross section is provided in Table 1. Beyond these two major families, a variety of other schemes have been applied including finite-element methods (e.g., Szymkiewicz, 1991; Venutelli, 2003), finite-volume methods that do not use the Godunov approach (e.g., Audusse et al., 2004, 2016; Catella et al., 2008; Katsaounis et al., 2004; Mohamed, 2014; Vazquez-Cendon, 1999; Xing and Shu, 2011; Ying et al., 2004), and finite-difference methods that do not apply the Preissmann scheme (e.g., Abbott and Ionescu, 1967; Aricò and Tucciarelli, 2007; Buntina and Ostapenko, 2008; Schippa and Pavan, 2008; Tucciarelli, 2003; Wang et al., 2000). A recent development is the introduction of discontinuous Galerkin (DG) methods, which can be thought of as a higher-order Godunov method (e.g., Kesserwani et al., 2008, 2009; Lai and Khan, 2012; Xing, 2014; Xing and Zhang, 2013).

It is clear from the above that there is no consensus on the best method for solving the SVEs. For high-resolution modeling that correctly preserves shocks and transcritical flows, it can be reasonably argued that finite-volume and DG schemes are more successful than finite-difference schemes. Beyond that broad observation, the question of whether a finite-volume method with a Godunov-like approach is better than a non-Godunov approach does not have a clear answer either in terms of accuracy or computational run times. However, in terms of large-scale systems the Preissmann scheme and finite-difference methods presently reign supreme with the ability to solve more than 15×10³ km of river on a desktop computer (Liu and Hodges, 2014). Nevertheless, we can see that a conservative finite-volume approach for large-scale systems would be preferred as the basis for simulating continental river dynamics (Hodges, 2013) as well as for the challenges of urban drainage modeling for the megacities that are growing across the earth. Being able to control numerical dissipation of momentum and ensure conservative fluxes will be increasingly important as advancing computing power pushes down the practical model grid scales.

Canelon (2009)Bollermann et al. (2013)Casas et al. (2010)Delis and Skeels (1998)Castellarin et al. (2009)Delis et al. (2000b)Chau and Lee (1991)Delis et al. (2000a)Chen et al. (2005)Ferreira et al. (2012)Gasiorowski (2013)Gottardi and Venutelli (2003)Islam et al. (2008)Goutal and Maurel (2002)Leandro and Martins (2016)Greenberg and Leroux (1996)Liu and Hodges (2014)Guinot (2009)Lyn and Altinakar (2002)Ivanova et al. (2017)Paiva et al. (2011, 2013)Kesserwani et al. (2010)Paz et al. (2010)Kurganov and Petrova (2007)Rosatti et al. (2011)Li and Chen (2006)Saavedra et al. (2003)Liang and Marche (2009)Sart et al. (2010)Monthe et al. (1999)Saleh et al. (2013)Pu et al. (2012)Sen and Garg (2002)Sanders (2001)Venutelli (2002)Sanders et al. (2003)Wu et al. (2004)Venutelli (2006)Zeng and Beck (2003)Wu and Wang (2007)Zhu et al. (2011)Ying and Wang (2008)

Table 1 Two decades of Preissman vs. Godunov.

3.3 The well-balanced problem and S₀

That finite-volume solutions are not commonly used in large-scale hydrologic and urban drainage models is a testament to their complexity. The difficulties associated with finite-volume solutions using the Cunge–Liggett conservative form of the SVEs have engendered a broad literature on well-balanced schemes – also known as schemes maintaining the “𝒞-property” – as derived in the study of Greenberg and Leroux (1996). A principal feature of a well-balanced scheme is that it provides exactly steady solutions for exactly steady flows. The most trivial requirement (which is often not met by unbalanced schemes) is that a flat free surface should result in exactly zero velocities. This problem is readily illustrated by considering Eq. (4) with Q=0, which implies S_f=0. Achieving the simple result of a flat free surface for Q=0 with the Cunge–Liggett form requires

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle \frac{\partial I_{\mathrm{1}}}{\partial x}}=I_{\mathrm{2}}+A{S}_{\mathrm{0}}\iff H+z_{\mathrm{b}}=\mathrm{constant},\end{array}$$

which implies the Cunge–Liggett form is only well-balanced if the geometry meets the following identity at every possible water surface level:

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial x}}& {\displaystyle }\underset{z_{\mathrm{b}}\left(x\right)}{\overset{z_R}{\int }}\left(z_R-z\right)B\mathrm{d}z-\underset{z_{\mathrm{b}}\left(x\right)}{\overset{z_R}{\int }}\left(z_R-z\right){\displaystyle \frac{\partial B}{\partial x}}\mathrm{d}z\\ \text{(13)}& {\displaystyle }& {\displaystyle }=A{S}_{\mathrm{0}}:z_{\mathrm{b}}<z_R\le {\mathit{\eta }}_{\max },\end{array}$$

where η_max is the maximum water surface elevation, z_b(x) is the local channel bottom elevation, and z_R encompasses all possible water surface elevations. Clearly, designing a numerical scheme that exactly preserves this relationship for nonuniform channels is a challenge, as evidenced by the breadth and complexities of studies focused on this issue (e.g., Audusse et al., 2004; Bollermann et al., 2013; Bouchut and Morales de Luna, 2010; Castro Diaz et al., 2007; Crnković et al., 2009; Kesserwani et al., 2010; Kurganov and Petrova, 2007; Li et al., 2017; Liang and Marche, 2009; Perthame and Simeoni, 2001; Xing, 2014). Failure to satisfy the well-balanced criteria results in models that generate spurious velocities; i.e., a mismatch in Eq. (13) indicates that the numerics provide momentum sources and sinks that are functions of channel shape and discretization rather than flow physics. An interesting approach to this problem was developed by Schippa and Pavan (2008), where Eq. (12) is used to replace I₂+AS₀ in the source term with $\partial I_{\mathrm{1}}/\partial x$ evaluated for a horizontal surface. Their approach ensures that any discretization will be well-balanced for a zero-velocity flow.

The work of Schippa and Pavan (2008) and review of other works on well-balanced schemes provide us with a key insight: the principal challenge for obtaining a well-balanced method is the channel bottom slope, S₀, which is often sharply varying or even discontinuous in a natural system. Furthermore, as a geometrical property, S₀ should be independent of the cross-sectional flow area (A), and yet is forced to be discretely related through Eq. (12). If we take this idea a step further, we can argue that the fundamental problem with the Cunge–Liggett form is that the physical forces that alter momentum (gravitational potential and hydrostatic pressure) are arbitrarily separated so that one is wholly within the source term and the other has an ad hoc split between conservative flux and source terms. Thus, we return to the idea put forward in the introduction that we should consider the free-surface elevation (piezometric pressure) instead of the water depth (hydrostatic pressure) as our primary forcing gradient.

In the next section, we shall see how the idea of shifting portions of the total piezometric pressure from source to flux can be used to develop a rigorous, conservative, and well-balanced finite-volume form of the SVEs that is simpler than those based on the Cunge–Liggett form.

4 Finite-volume SVEs with minimal approximations

4.1 Continuity

Although we are focused on the momentum equation, for completeness we will start with continuity. The general arrangement of the control volume for an irregular channel and the vectors used in the following discussion are illustrated in Fig. 1. Applying only the incompressibility approximation for a uniform-density fluid, the volume-integrated continuity equation is

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle \frac{\partial V}{\partial t}}+\underset{S}{\oint }{\mathit{u}}_k{n}_k\mathrm{d}A=S_{\mathrm{V}},\end{array}$$

where the Einstein summation convention is applied on repeated subscripts, u_k is a vector velocity, n_k is a unit normal vector defined as positive pointing outward from a control volume V, and S_V is a volume source (S_V<0 for a sink).

Figure 1 General arrangement of control-volume element (V_e) and its neighbors for the irregular channel. Unit normal vectors n_k are always perpendicular to cross-sectional areas (A_u, A_d) and pointing outward from control volume. The element length (L_e) is measured along the channel. Unit vectors ${\widehat{\mathit{i}}}_k$ are coincident with the free-surface slope in the streamwise direction and can be defined as local continuous functions. The velocity vector u_k is approximated as parallel to ${\widehat{\mathit{i}}}_k$. The angle measured from n_k to ${\widehat{\mathit{i}}}_k$ is ψ. The free-surface elevations η_u, η_e, and η_d are cross-section uniform elevations at the upstream face, for the element center, and the downstream face, respectively. Note that volume is V throughout the derivations with u and U used for continuous and spatially averaged velocities, respectively.

A semi-discrete finite-volume representation of continuity can be directly written as

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle \frac{\partial V_{\mathrm{e}}}{\partial t}}=Q_{\mathrm{u}}-Q_{\mathrm{d}}+q_{\mathrm{e}}{L}_{\mathrm{e}},\end{array}$$

where Q and L represent the flow rate and element length, and subscripts “e”, “u”, and “d” denote characteristic values for the control-volume element, nominal upstream face, and nominal downstream face, respectively. Here we use the nominal flow direction as the global downhill direction of the channel that is assigned at the network level. The average lateral inflow per unit length is q_e, and a flow Q>0 is from the upstream to downstream direction. Reversals of flow from the nominal flow direction are handled with Q<0.

4.2 Momentum

The control-volume form of the Navier–Stokes momentum equation in a direction defined by unit vector $\widehat{\mathit{i}}$ in a Cartesian frame is

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial t}}\underset{V}{\int }{u}_{\widehat{\mathit{i}}}\mathrm{d}V& {\displaystyle }+\underset{S}{\oint }{u}_{\widehat{\mathit{i}}}{u}_k{n}_k\mathrm{d}A+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{S}{\oint }p{\widehat{i}}_k{n}_k\mathrm{d}A\\ \text{(16)}& {\displaystyle }& {\displaystyle }=\underset{S}{\oint }\mathit{\nu }{\displaystyle \frac{\partial u_j}{\partial x_k}}{\widehat{i}}_j{n}_k\mathrm{d}A+\underset{V}{\int }{g}_{\widehat{\mathit{i}}}\mathrm{d}V,\end{array}$$

where $u_{\widehat{\mathit{i}}}$ is a velocity vector component in the $\widehat{\mathit{i}}$ direction (i.e., the direction that is a priori of interest), u_k represents velocity components along Cartesian axes, the component of the gravity vector in the $\widehat{\mathit{i}}$ direction is $g_{\widehat{\mathit{i}}}$, the kinematic viscosity is ν, and p represents the thermodynamic pressure. Note that this formulation can be related to any arbitrary Cartesian axes. In many common derivations, $\widehat{\mathit{i}}$ is approximated as coincident with an x axis that is a horizontal vector in the streamwise direction. In the following, we will show that this approximation is not required. Instead, we treat this as a simplification that can be applied to the final equation form.

4.3 Advection terms

The $\widehat{\mathit{i}}$ direction for momentum, Eq. (16), is a vector associated with the $u_{\widehat{\mathit{i}}}$ velocity, which is not necessarily coincident with the normal vector n_k at a flux surface of a finite volume (in contrast to the case where $\widehat{\mathit{i}}$ is taken as horizontal). For a gradually varying open-channel flow we can take the $\widehat{\mathit{i}}$ vector as the nominal downstream direction along the channel centerline described by a vector that lies along the free surface, as illustrated in Fig. 1. Thus, this vector is local (as opposed to being forced into coincidence with a Cartesian axis) and changes along the channel with the slope of the free surface. It follows that a discrete control-volume formulation developed from Eq. (16) can be globally exact as V_e→0. In contrast, a derivation that takes $\widehat{\mathit{i}}$ as a vector in the horizontal direction has a momentum conservation error proportional to cos ψ, where ψ is an angle between the horizontal vector and the free surface; such an error does not vanish as V_e→0 unless the free surface is flat across the length of the element. This idea helps illustrate one of the subtle implications of the Godunov approach in which the channel is imagined as having a piecewise flat free surface: cos ψ=1 is then an identity within the approximation of the physics rather than an approximation within the mathematics.

For Eq. (16), the upstream element face is required to be vertical and normal to the smooth channel centerline in a horizontal plane, as illustrated in Fig. 1. The free surface at the centerline has an angle of ψ(x) to the horizontal so that the discrete nonlinear momentum term in the $\widehat{\mathit{i}}$ direction across the upstream face is formally

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }\underset{A_{\mathrm{u}}}{\int }{u}_{\widehat{\mathit{i}}}{u}_k{n}_k\mathrm{d}A={\left(\mathit{\beta }Q{U}_{\widehat{\mathit{i}}}\right)}_{\mathrm{u}},\end{array}$$

where Q is the flow rate, $U_{\widehat{\mathit{i}}}$ is the average streamwise velocity over the cross-sectional area, the “u” subscript indicates the upstream face (rather than vector components), and β is the momentum coefficient for the streamwise velocity, defined as

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }\mathit{\beta }\equiv {\displaystyle \frac{\mathrm{1}}{A{\left(U_{\widehat{\mathit{i}}}\right)}^{\mathrm{2}}}}\underset{A}{\int }{\left(u_{\widehat{\mathit{i}}}\right)}^{\mathrm{2}}\mathrm{d}A.\end{array}$$

Note that the only approximation in the convective term of Eq. (17) is that the streamwise velocity is parallel to the free surface. However, the interpretation of the $u_{\widehat{\mathit{i}}}$ and u_kn_k terms may not be obvious, so further explanation is provided in Appendix A. For notational convenience, it is useful to let $U\equiv U_{\widehat{\mathit{i}}}$. However, since $Q=A\int u_k{n}_k\mathrm{d}A$, strictly speaking this requires an unconventional Q=AUcos ψ. If the channel is straight and the free surface is linear so that the upstream face is parallel to the downstream face and there is a single value of ψ, it follows that an exact finite-volume integration of the nonlinear advection term is

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }\underset{S}{\oint }{u}_{\widehat{\mathit{i}}}{u}_k{n}_k\mathrm{d}A=-{\mathit{\beta }}_{\mathrm{u}}{Q}_{\mathrm{u}}{U}_{\mathrm{u}}+{\mathit{\beta }}_{\mathrm{d}}{Q}_{\mathrm{d}}{U}_{\mathrm{d}}-M_{\mathrm{e}},\end{array}$$

where M_e represents any integrated sources (+) and sinks (−) of momentum per unit mass in the finite-volume element associated with the q_eL_e lateral fluxes of Eq. (15). Note that M_e terms are typically neglected in SVE solvers. For the more general case where the channel is curved between the upstream and downstream faces and the free-surface gradient changes (as in Fig. 1), the above integration becomes an approximation that is only exactly satisfied in the limit as the control-volume length goes to zero. For present purposes, the use of a gradually varying streamwise $\widehat{\mathit{i}}$ direction implies that pressure is perfectly redirecting momentum through bends and aligning the momentum with the free surface. These are (generally) unstated approximations used in common 1-D SVE formulations. However, it should be noted that this perfect momentum redirection is not precisely correct; e.g., secondary circulation in bends affects bed shear, velocity distribution, and frictional losses (Blanckaert and Graf, 2004). Arguably, losses associated with flow redirection in channel bends and/or rapid changes in the free-surface gradient should be built into the frictional term in any model. Unfortunately, this remains a relatively poorly studied area at the interface of hydrology and hydraulics. The curvature effects on the equations can be written as perturbation terms that relate the channel width to the radius of curvature (Hodges and Imberger, 2001; Hodges and Liu, 2014), but these ideas have yet to be exploited in developing curvature effects in SVE numerical models.

4.4 Pressure decomposition

For the pressure term in Eq. (16) we follow the traditional approach for incompressible flows of defining a modified pressure ($\stackrel{\mathrm{̃}}{P}$) that includes the gravity term, which requires $\stackrel{\mathrm{̃}}{P}=p+\mathit{\rho }gz$. More formally we define

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{S}{\oint }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A\equiv {\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{S}{\oint }p{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A-\underset{V}{\int }{g}_{\widehat{\mathit{i}}}\mathrm{d}V.\end{array}$$

The surface integral for the modified pressure decomposes into integrations over the upstream cross section (A_u), the downstream cross section (A_d), the channel bottom (A_B), and the free surface (A_η) as

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{S}{\oint }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\left\{\underset{A_{\mathrm{u}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A+\underset{A_{\mathrm{d}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A\right.\\ \text{(21)}& {\displaystyle }& {\displaystyle }\left.+\underset{A_{\mathrm{B}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A+\underset{A_{\mathit{\eta }}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A\right\}.\end{array}$$

Note that A_B includes both the bottom and side walls (if any). The above is an exact pressure term without imagining the geometry to be rectilinear or that the free surface has any simplified shape. The only geometric requirement of the volume is that faces A_u and A_d must be vertical planes that cut across the channel, which is necessary for consistency with the convective term. The first simplification we can introduce is to approximate the free surface as uniform at any cross section. The slope of the free surface is then coincident with the $\widehat{\mathit{i}}\left(x\right)$ vector everywhere, which is aligned with the free surface at the channel centerline. With the free surface aligned with $\widehat{\mathit{i}}\left(x\right)$, the modified pressure is normal to the free surface everywhere and cannot contribute to the streamwise momentum (which we have defined as parallel to $\widehat{\mathit{i}}\left(x\right)$ in deriving the advection terms, above). It follows that $\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k$ is identically zero at the free surface and the last term in Eq. (21) vanishes.

4.5 Piezometric and non-hydrostatic pressure

It is convenient to introduce a decomposition of the modified pressure ($\stackrel{\mathrm{̃}}{P}$) into a piezometric pressure (P) and non-hydrostatic pressure ($\stackrel{\mathrm{˘}}{P}$), where only the latter is nonuniform over a cross section. The non-hydrostatic pressure is defined by the difference between the modified pressure and the piezometric pressure:

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{˘}}{P}(x,y,z)\equiv \stackrel{\mathrm{̃}}{P}(x,y,z)-P\left(x\right).\end{array}$$

The piezometric pressure is formally the sum of the hydrostatic pressure and the gravitational potential at any point z for $z_{\mathrm{b}}\le z\le \mathit{\eta }$, which provides P(x, y, z)≡ρg[η(x, $y)-z]+\mathit{\rho }gz=\mathit{\rho }g\mathit{\eta }(x$, y). For the SVEs, the free-surface elevation can be considered uniform over the channel breadth (i.e., neglecting cross-channel tilt in channel bends). It follows that the piezometric pressure is

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }P\left(x\right)=\mathit{\rho }g\mathit{\eta }\left(x\right),\end{array}$$

which is uniform over a vertical cross section. The non-hydrostatic pressure was neglected by de Saint-Venant and arguably should be neglected in any 1-D momentum equation bearing his name. However, for completeness we retain the non-hydrostatic pressure terms in the derivation below but neglect them in discrete forms of the piezometric pressure terms in Sect. 5.

4.6 Pressure on flow faces

To further simplify the first two pressure terms in Eq. (21), we define ψ(x) as the angle between a horizontal line and the free surface, η(x), measured clockwise positive from the horizontal line pointing downstream. Thus, the $\widehat{\mathit{i}}$ vector is generally at an angle ±ψ(x) and pointing in the nominal downstream direction, as shown in Fig. 1. At the upstream cross section the pressure force without approximations is

$$\begin{array}{}\text{(24)}& {\displaystyle }{\displaystyle }\underset{A_{\mathrm{u}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A=-\left(\cos {\mathit{\psi }}_{\mathrm{u}}\right)\underset{A_{\mathrm{u}}}{\int }\stackrel{\mathrm{̃}}{P}(x,y,z)\mathrm{d}A,\end{array}$$

where subscript “u” indicates values at the upstream cross section. Similarly, the downstream cross section provides

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }\underset{A_{\mathrm{d}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A=+\left(\cos {\mathit{\psi }}_{\mathrm{d}}\right)\underset{A_{\mathrm{d}}}{\int }\stackrel{\mathrm{̃}}{P}(x,y,z)\mathrm{d}A,\end{array}$$

where subscript “d” indicates values at the downstream cross section. Note that, in the above and in the following derivations, the ubiquitous 1∕ρ coefficient of all pressure terms are omitted for clarity. Introducing the piezometric and non-hydrostatic pressure split of Eq. (22) provides

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }\underset{A_{\mathrm{u}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A=-P_{\mathrm{u}}{A}_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}-\left(\cos {\mathit{\psi }}_{\mathrm{u}}\right)\underset{A_{\mathrm{u}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A,\text{(27)}& {\displaystyle }& {\displaystyle }\underset{A_{\mathrm{d}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A=+P_{\mathrm{d}}{A}_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}+\left(\cos {\mathit{\psi }}_{\mathrm{d}}\right)\underset{A_{\mathrm{d}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A.\end{array}$$

Note that using the piezometric pressure instead of the hydrostatic pressure ensures that the only term requiring discrete integration at the upstream and downstream surfaces is the non-hydrostatic pressure. That is, PAcos ψ is not an approximate integral whose adequacy depends on simplifying assumptions in geometry (as is the case for integrals of hydrostatic pressure in the Cunge–Liggett form) but is instead an exact integration of piezometric pressure for any cross-section shape.

Figure 2 Pressure decomposition to obtain streamwise contribution.

4.7 Pressure on bottom topography

The third pressure term in Eq. (21) is more challenging than the pressure on the flow faces as the bottom surface normal (n_k) varies with irregular topography. Hence, the local piezometric pressure contribution in the streamwise direction ($\widehat{\mathit{i}}$) at any position (x, y) depends not only on the local water surface elevation but on the local irregularities in topography. It follows that integrating this term over a control volume requires some approximation of the subgrid topography. To arrive at a simpler formulation we note that the pressure acting normal to an arbitrary topography element at position (x, y) with surface normal n_k will have force components that can be resolved along a set of local Cartesian axes. One axis is taken along a topography slope angle of θ that lies in the same vertical plane as the streamwise direction $\widehat{\mathit{i}}$, as illustrated in Fig. 2. We can imagine Fig. 2 as an infinitesimal slice across a channel of irregular topography such that a series of these slices (with different θ) can represent any cross-section shape. The second local Cartesian axis is taken within the same vertical plane but perpendicular to the axis defined by θ. The third axis, perpendicular to the other two, is necessarily horizontal and in the cross-channel direction (out of the page in Fig. 2). We are interested only in how the topography contributes to pressure in the streamwise $\widehat{\mathit{i}}$ direction, so by definition the cross-channel pressure component is irrelevant. This approach is entirely consistent with changes in depth across the channel and changes in breadth along the channel – both merely alter the surface normal vector n_k that is resolved into the local Cartesian system based on the local slope direction of θ.

Figure 3 Stair-step approximation for pressure along the bottom where tread is parallel to the free surface and riser is perpendicular for linear slopes of both the bottom and free surface.

To isolate the pressure forces acting in the $\widehat{\mathit{i}}$ streamwise direction, as a conceptual model we can imagine the topography in the infinitesimal slice of Fig. 2 replaced with a set of m=1 … N stair steps, where the treads are locally parallel to the free surface and the risers are normal to the free surface, as illustrated in Fig. 3. Clearly, as N→∞ we will recover a continuous approximation so there is no need to actually consider the discrete stair steps in a solution method – the steps are merely to illustrate what is otherwise a mathematical abstraction in vector calculus. As illustrated in Fig. 4, we can imagine the stair-step risers as thin planar strips across the entire wetted perimeter that provide a discrete representation of the irregular channel cross-section structure. The only requirements for this conceptual model are that the free surface at longitudinal position x is uniform over the cross section and the slope is aligned with $\widehat{\mathit{i}}\left(x\right)$. This conceptual model could also be envisioned in 3-D as discrete rectilinear bricks that approximate the topography: the upper surface of the brick is always aligned parallel with the free surface, the side face is across the channel, and only the front (or back) face is perpendicular to $\widehat{\mathit{i}}$ and contributes a topographic pressure force in the streamwise direction. As the brick dimensions go to zero the continuous topography is obtained.

Figure 4 Three-dimensional conceptual model of stair-step riser (gray) as a 2-D planar area of A_R(m) separating two tread areas (blue, red).

Figure 5 Conceptual model of piecewise linear approximation of nonlinear free surface and stair-step approximation of nonlinear topography with the cross-sectional area monotonically increasing. Pressure on the discrete riser area A_R(m) is locally aligned with the free-surface slope directly above. The vertical scale and thus the tilt of A_R(m) are exaggerated for illustrative purposes.

Since the stair-step treads are (by definition) normal to the modified pressure above, it follows that the only pressure contributions to the momentum in the streamwise $\widehat{\mathit{i}}$ direction are on the risers, with individual areas A_R(m) for m=1 … N stair steps. Because the pressure contribution for increasing cross-sectional area (i.e., steps down as in Fig. 3) will be opposite of the pressure contribution for decreasing cross-sectional area (i.e., steps upward), it is convenient to introduce a function ${\mathit{\gamma }}_{\left(m\right)}=\pm \mathrm{1}$ to account for the change in sign needed for the direction of the pressure force. We can formally define ${\mathit{\gamma }}_{\left(m\right)}\equiv {\widehat{\mathit{i}}}_k{\widehat{j}}_k$ at step m where ${\widehat{j}}_k$ is the normal unit vector (pointing outwards) from the A_R(m) riser, as shown in Figs. 5 and 6 for two different nonlinear water surface profiles over identical bottom topography. It follows that

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }\underset{A_{\mathrm{B}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A\approx -\sum _{m=\mathrm{1}}^N{\mathit{\gamma }}_{\left(m\right)}\underset{A_{\mathrm{R}\left(m\right)}}{\int }\stackrel{\mathrm{̃}}{P}(x,y,z)\mathrm{d}A.\end{array}$$

The above summation can be written as an integral over the length L of the finite-volume element as N→∞.

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }\underset{A_{\mathrm{B}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A=-\underset{L}{\int }\mathit{\gamma }\left(x\right)\underset{A_{\mathrm{R}}\left(x\right)}{\int }\stackrel{\mathrm{̃}}{P}(x,y,z)\mathrm{d}A\mathrm{d}x,\end{array}$$

where $\mathit{\gamma }\left(x\right)={\widehat{\mathit{i}}}_k{\widehat{j}}_k$ is the continuous counterpart to the discrete γ_(m). Note that this conceptual model is valid even for non-monotonic behavior of the riser area (e.g., Fig. 6) as long as the A_R(m)'s are continuous and smooth as N→∞. However as discussed in Sect. 5, below, extremes of non-monotonic behavior can make it difficult to create a consistent discrete equation for the topographic pressure for a control volume of finite size.

Figure 6 Conceptual model of piecewise linear approximation of nonlinear free surface and stair-step approximation of nonlinear bottom with cross-sectional area increasing from A_u to the center and then decreasing from the center to A_d. Pressure on the discrete riser area A_R(m) is locally aligned with the free-surface slope directly above. The vertical scale and thus the tilt of A_R(m) are exaggerated for illustrative purposes. Note that the topography is identical to Fig. 5 but the stair steps are different due to the alteration of the free surface.

For further simplification, it is convenient to introduce the piezometric/non-hydrostatic splitting of Eq. (22) where, by definition, the piezometric pressure is uniform over a vertical cross section (e.g., a stair-step riser). Let P_(m) represent the piezometric pressure at the m stair-step riser with area A_R(m) so that over a control volume for steps m=1 … N. It follows that

$$\begin{array}{ll}{\displaystyle }\sum _{m=\mathrm{1}}^N\underset{A_{\mathrm{R}\left(m\right)}}{\int }{P}_{\left(m\right)}\mathrm{d}A& {\displaystyle }=\sum _{m=\mathrm{1}}^N{P}_{\left(m\right)}\underset{A_{\mathrm{R}\left(m\right)}}{\int }\mathrm{d}A\\ \text{(30)}& {\displaystyle }& {\displaystyle }=\underset{L}{\int }P\left(x\right)A_{\mathrm{R}}\left(x\right)\mathrm{d}x\mathrm{for}N\to \mathrm{\infty },\end{array}$$

where the continuous A_R(x) represents the effective bottom area contributing to the piezometric pressure force in the streamwise $\widehat{\mathit{i}}\left(x\right)$ direction. Applying this continuous A_R(x) form and the hydrostatic/non-hydrostatic splitting to Eq. (29), we obtain

$$\begin{array}{ll}{\displaystyle }\underset{A_{\mathrm{B}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A& {\displaystyle }=-\underset{L}{\int }\mathit{\gamma }\left(x\right)P\left(x\right)A_{\mathrm{R}}\left(x\right)\mathrm{d}x\\ \text{(31)}& {\displaystyle }& {\displaystyle }-\underset{L}{\int }\mathit{\gamma }\left(x\right)\underset{A_{\mathrm{R}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A\mathrm{d}x,\end{array}$$

which is a complete representation of the topographic pressure contribution to streamwise momentum in terms of the effective riser areas – i.e., the contribution based on the component of the bottom normal projected in the streamwise direction.

The stair-step conceptual model and A_R allow us to consider the pressure effects along the $\widehat{\mathit{i}}$ direction due to the changing cross-sectional area of the channel without introducing the separate force terms I₁ and I₂ of the Cunge–Liggett form of the SVEs. An interesting part of this model is that A_R(x) over a control volume is a function of both the local bottom topography and the local free-surface slope. That is, comparing Figs. 5 and 6 we see that A_R(m) riser areas are different, despite the identical bottom topography. Returning to our idea that 3-D topography can be represented by discrete bricks, we can imagine each brick is pinned on an axis that allows it to locally rotate to different angles so that the upper surface is always parallel to the free surface. Again, the continuous topography is recovered as the brick size goes to zero, but the brick rotation allows the representation of the different topography effects that are caused by the interaction between the change in relationship between the downstream vector, $\widehat{\mathit{i}}$, and the bottom normal vector as the free-surface profile evolves. Thus, A_R is a dynamic representation of the interaction between the free-surface gradient and bottom topography that controls the effective along-stream pressure gradient of converging or diverging flow areas.

In theory, we might directly compute ∫PA_Rdx over a control volume; however, it seems likely that direct discretization of subgrid topography could cause unbalanced momentum source terms. In effect, computing ∫PA_Rdx has the same complexity as computing the I₁ or I₂ terms of the Cunge–Liggett form and gains us little. Thus, it is useful to consider limiting approximations that can be developed from examining the geometry of Figs. 3 and 5. In the simplest case where both the free surface and bottom have linear slopes (Fig. 3), geometry for the infinitesimal slice requires that

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }{A}_{\mathrm{d}}\cos \mathit{\psi }-A_{\mathrm{u}}\cos \mathit{\psi }=\sum _{m=\mathrm{1}}^N{A}_{\mathrm{R}\left(m\right)},\end{array}$$

where A_dcos ψ represents the cross-sectional area normal to the free surface at the downstream cross section and there is only one value of ψ over the control volume for the linear free-surface slope. Note that the above is an identity for any discrete N stair steps for the linear case. The geometry is somewhat more complex for the nonlinear system of Fig. 5, but the m stair step with free-surface angle ψ_(m) must satisfy a similar geometric identity:

$$\begin{array}{}\text{(33)}& {\displaystyle }{\displaystyle }{A}_{(m-\mathrm{1}/\mathrm{2})}\cos {\mathit{\psi }}_{\left(m\right)}-A_{(m+\mathrm{1}/\mathrm{2})}\cos {\mathit{\psi }}_{\left(m\right)}=A_{\mathrm{R}\left(m\right)},\end{array}$$

where $A_{(m\pm \mathrm{1}/\mathrm{2})}\cos {\mathit{\psi }}_{\left(m\right)}$ represents the areas normal to the free surface on the upstream and downstream edges of the m piecewise linear stair step. For the nonlinear free surface and/or topography over adjacent linear stair steps there is a discontinuity of the treads for adjacent steps as $\cos {\mathit{\psi }}_{(m-\mathrm{1})}\ne \cos {\mathit{\psi }}_{\left(m\right)}$, so the discrete summation of the stair-step areas over a control volume provides an approximation rather than an identity:

$$\begin{array}{}\text{(34)}& {\displaystyle }{\displaystyle }{A}_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}\approx \sum _{m=\mathrm{1}}^N{A}_{\mathrm{R}\left(m\right)}.\end{array}$$

However, in the limit as N→∞ we have a single value of cos ψ at any point along a smooth free surface so that the continuous form provides an identity:

$$\begin{array}{}\text{(35)}& {\displaystyle }{\displaystyle }{A}_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}=\underset{L}{\int }{A}_{\mathrm{R}}\left(x\right)\mathrm{d}x.\end{array}$$

To generalize the above for A_d<A_u, we can use the $\mathit{\gamma }\left(x\right)=\pm \mathrm{1}$ that was introduced for the pressure direction in Eq. (29). Values of $\mathit{\gamma }\left(x\right)=+\mathrm{1}$ indicate that the cross-section area is increasing across location x in the streamwise direction (as in Figs. 3 and 5), whereas $\mathit{\gamma }\left(x\right)=-\mathrm{1}$ indicates that the cross-sectional area is decreasing (as in the latter portion of Fig. 6). It follows that

$$\begin{array}{}\text{(36)}& {\displaystyle }{\displaystyle }{A}_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}=\underset{L}{\int }\mathit{\gamma }\left(x\right)A_{\mathrm{R}}\left(x\right)\mathrm{d}x\end{array}$$

is an identity that should be satisfied for any control volume where the bottom topography and free surface are continuous and smooth. Note that if Fig. 3 is imagined as one of many infinitesimal slices (with varying θ) that make up a channel cross section, it should be obvious that Eq. (36) also applies for a finite volume with irregular (smooth) topography.

To handle the integration of A_R(x) in the piezometric pressure term of Eq. (31), we introduce a quadrature function λ(x), defined as

$$\begin{array}{}\text{(37)}& {\displaystyle }{\displaystyle }\mathit{\lambda }\left(x\right)\equiv {\displaystyle \frac{\mathit{\gamma }\left(x\right)A_{\mathrm{R}}\left(x\right)}{A_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}}}.\end{array}$$

Note that with Eq. (36), this implies the identity

$$\begin{array}{}\text{(38)}& {\displaystyle }{\displaystyle }\underset{L}{\int }\mathit{\lambda }\left(x\right)\mathrm{d}x=\mathrm{1}.\end{array}$$

Using the above in the first term on the RHS of Eq. (31), we obtain

$$\begin{array}{ll}{\displaystyle }\underset{L}{\int }\mathit{\gamma }\left(x\right)P\left(x\right)A_{\mathrm{R}}\left(x\right)\mathrm{d}x=& {\displaystyle }\left(A_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}\right)\\ \text{(39)}& {\displaystyle }& {\displaystyle }\times \underset{L}{\int }P\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x.\end{array}$$

Thus, the introduction of λ allows us to extract a multiplier from the control-volume integral of the bottom pressure. As a result, λ(x) is merely a distribution, or weighting function, for integration of P(x). The full bottom pressure term, Eq. (31), can be written as

$$\begin{array}{ll}{\displaystyle }\underset{A_{\mathrm{B}}}{\int }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A& {\displaystyle }=-\left(A_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}\right)\underset{L}{\int }P\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x\\ \text{(40)}& {\displaystyle }& {\displaystyle }-\underset{L}{\int }\mathit{\gamma }\left(x\right)\underset{A_{\mathrm{R}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A\mathrm{d}x.\end{array}$$

Note that λ weighting cannot be readily applied to the non-hydrostatic term because the non-hydrostatic pressure on the bottom has spatial distributions in both the vertical and cross-channel directions that cannot be assumed to be negligible; hence we cannot pass $\stackrel{\mathrm{˘}}{P}$ through the A_R integration as was done in Eq. (30) for P.

We can think of λ(x) as a weighting function of the conceptual stair-step riser areas over the control-volume length, which controls where the piezometric pressure gradients have their greatest effect. For example, in Fig. 3 the stair-step risers are uniformly distributed such that we can use $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$, which meets the identity requirement of Eq. (38). In contrast, Fig. 5 implies that λ(x) is perhaps a quadratic function. Figure 6 presents a challenge as λ(x) should reverse in sign between the upstream and downstream faces. A key point in the new finite-volume derivation is that λ(x) controls the representation of the free surface and bottom topography within a volume. This can be contrasted to the Godunov approach that uses piecewise linear approximations for both the bottom and free-surface elevations (see Sect. 3). Several discrete approaches to the approximation of λ(x) are examined in Sect. 5, although the full consequences and utility of the λ approach will require more extensive investigation for both theoretical limitations and practical discretization schemes.

4.8 Combining pressure terms

In summary, the pressure terms of Eq. (21) can be written using Eqs. (26), (27), and (40), resulting in

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{S}{\oint }\stackrel{\mathrm{̃}}{P}{\widehat{\mathit{i}}}_k{n}_k\mathrm{d}A& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{A}_{\mathrm{u}}{P}_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{A}_{\mathrm{d}}{P}_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\left(A_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}\right)\underset{L}{\int }P\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\cos {\mathit{\psi }}_{\mathrm{u}}}{\mathit{\rho }}}\underset{A_{\mathrm{u}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\cos {\mathit{\psi }}_{\mathrm{d}}}{\mathit{\rho }}}\underset{A_{\mathrm{d}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A\\ \text{(41)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{L}{\int }\mathit{\gamma }\left(x\right)\underset{A_{\mathrm{R}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A\mathrm{d}x,\end{array}$$

where the last three terms are the non-hydrostatic pressure effects that are typically neglected in the SVEs.

4.9 Viscous term

The remaining term in Eq. (16) is the viscous term, which is treated as an empirical function in all but the most highly resolved models of simple systems – note that Decoene et al. (2009) provide a comprehensive and rigorous approach for friction that has not yet been fully considered in SVE models. For the present purposes, we will retain the simple friction slope form with an assumption of uniform behavior over space, i.e.,

$$\begin{array}{}\text{(42)}& {\displaystyle }{\displaystyle }\underset{S}{\oint }\mathit{\nu }{\displaystyle \frac{\partial u_j}{\partial x_k}}{\widehat{\mathit{i}}}_j{n}_k\mathrm{d}A=-g\underset{V_{\mathrm{e}}}{\int }{S}_{\mathrm{f}}\left(x\right)\mathrm{d}V\approx -g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)},\end{array}$$

where S_f(e) is the average friction slope over the control volume V_e.

4.10 Finite volume for momentum

Putting together the above, Eq. (16) can be written in a finite-volume form as

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial t}}\left(U_{\mathrm{e}}{V}_{\mathrm{e}}\right)& {\displaystyle }={\mathit{\beta }}_{\mathrm{u}}{Q}_{\mathrm{u}}{U}_{\mathrm{u}}-{\mathit{\beta }}_{\mathrm{d}}{Q}_{\mathrm{d}}{U}_{\mathrm{d}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{A}_{\mathrm{u}}{P}_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{A}_{\mathrm{d}}{P}_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\left(A_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}\right)\\ {\displaystyle }& {\displaystyle }\times \underset{L}{\int }P\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x-{\displaystyle \frac{\cos {\mathit{\psi }}_{\mathrm{u}}}{\mathit{\rho }}}\underset{A_{\mathrm{u}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\cos {\mathit{\psi }}_{\mathrm{d}}}{\mathit{\rho }}}\underset{A_{\mathrm{d}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{L}{\int }\mathit{\gamma }\left(x\right)\\ \text{(43)}& {\displaystyle }& {\displaystyle }\underset{A_{\mathrm{R}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A\mathrm{d}x-g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)}+M_{\mathrm{e}},\end{array}$$

where U_e is the element velocity in the streamwise direction, V_e is the element volume, and the relationship between U and Q is given by

$$\begin{array}{}\text{(44)}& {\displaystyle }{\displaystyle }Q=AU\cos \mathit{\psi }.\end{array}$$

Note that Q>0 and U>0 imply flow in the nominal downstream direction, whereas Q<0 and U<0 imply flow in the nominal upstream direction. At this point we have introduced only four approximations: (1) uniform-density incompressibility, (2) the effect of channel curvature is either negligible or handled in an empirical viscous term, (3) the cross-channel variability in the free-surface slope is negligible, and (4) a friction-slope model can be used to represent integrated viscous effects over a control volume. In addition we have a geometric restriction that the upstream and downstream control-volume cross sections must be vertical planes that are orthogonal (in the horizontal plane) to the mean flow direction.

For convenience in exposition, for the remainder of this paper we will apply the hydrostatic approximation ($\stackrel{\mathrm{˘}}{P}=\mathrm{0}$) along with approximations for the small slope (cos ψ≈1) and uniform cross-section velocity (β≈1). Furthermore, we will limit our focus to flows without lateral momentum sources (M_e=0). These simplifications allow us to focus attention on the pressure source term, which is the primary new contribution in this derivation. The resulting simplification of Eq. (43) can be presented with conservative terms on the LHS and source terms on the RHS as

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial t}}\left(U_{\mathrm{e}}{V}_{\mathrm{e}}\right)& {\displaystyle }-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{d}}}}-g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}+g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}\\ \text{(45)}& {\displaystyle }& {\displaystyle }=g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\underset{L}{\int }\mathit{\eta }\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x-g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)},\end{array}$$

where definition of piezometric pressure, Eq. (23), is used to substitute P=ρgη. A more formal finite-volume integral presentation would be

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial t}}\underset{V}{\int }u\mathrm{d}V& {\displaystyle }-\underset{A_{\mathrm{u}}}{\int }{u}^{\mathrm{2}}\mathrm{d}A+\underset{A_{\mathrm{d}}}{\int }{u}^{\mathrm{2}}\mathrm{d}A-g\underset{A_{\mathrm{u}}}{\int }\mathit{\eta }\mathrm{d}A+g\underset{A_{\mathrm{d}}}{\int }\mathit{\eta }\mathrm{d}A\\ \text{(46)}& {\displaystyle }& {\displaystyle }=g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\underset{L}{\int }\mathit{\eta }\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x-g\underset{V}{\int }{S}_{\mathrm{f}}\mathrm{d}V.\end{array}$$

Equation (45) can be reduced to a differential equation using V=AL as L→0. Dividing through by L and taking the limit as L→0 provides

$$\begin{array}{}\text{(47)}& {\displaystyle }{\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{Q^{\mathrm{2}}}{A}}+gA\mathit{\eta }\right)=g\mathit{\eta }{\displaystyle \frac{\partial A}{\partial x}}-gA{S}_{\mathrm{f}}.\end{array}$$

If we substitute geometric identities $\mathit{\eta }\equiv H+z_{\mathrm{b}}$ and $S_{\mathrm{0}}\equiv -\partial z_{\mathrm{b}}/\partial x$, we see that the above becomes identical to Eq. (8). Thus, our finite-volume derivation is exactly consistent with the commonly used differential SVEs that is posed using S₀. However, from Eqs. (45) and (47) we see that the free-surface source term in the differential form, $g\mathit{\eta }\partial A/\partial x$, is related to the more interesting integral source term in the finite-volume form:

$${\displaystyle }{\displaystyle }g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\underset{L}{\int }\mathit{\eta }\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x.$$

This piezometric pressure term can be thought of as the integrated free-surface/topography effects over a control volume of finite size. It is clear that this term collapses to $g\mathit{\eta }\partial A/\partial x$ as L→0, but the integral form cannot be readily inferred from the differential form. Through approximations of this integral term we can obtain a variety of different finite-volume forms of the SVEs, as outlined in the following section.

5 Approximate finite-volume forms of the SVEs

5.1 General form

We are interested in approximate forms of the finite-volume SVEs that arise from discretization choices in the integral pressure source term derived above, so for exposition it is convenient to start from the approximate form in Eq. (45) with V_e=A_eL_e and Q_e=A_eU_e. These approximations can be readily reversed to provide a more complete equation, but the simple form for further analysis is

$$\begin{array}{}\text{(48)}& {\displaystyle }{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial Q_{\mathrm{e}}}{\partial t}}-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{u}}}}-g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}+g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}=T_{\mathrm{e}}-g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)},\end{array}$$

where T_e is the source term for integrated free-surface/topography pressure effects obtained from the RHS of Eq. (45):

$$\begin{array}{}\text{(49)}& {\displaystyle }{\displaystyle }{T}_{\mathrm{e}}\equiv g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\underset{L}{\int }\mathit{\eta }\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x.\end{array}$$

Note that the LHS of Eq. (48) is discretely conservative in that a summation over all elements will cause all the LHS terms to identically vanish except for the $\partial /\partial t$ and boundary conditions. Thus, the RHS shows the source terms; i.e., the traditional friction term and a term representing nonlinearities in the free surface interacting with topography can create and destroy momentum. The system is inherently well-balanced as discussed in Sect. 3, as long as T_e identically balances the spatially integrated piezometric pressure gradient gA_dη_d−gA_uη_u for a flat free surface.

Equation (49) admits a wide variety of approximate equations, depending on the form chosen for the quadrature of η(x) and λ(x) over the length of a finite volume. Arguably, a simple finite-volume discrete method will have three values of η that characterize a control volume: η_u, η_e, and η_d, as illustrated in Fig. 1. Different numerical methods can be constructed by using different approximations constructed from these three values. Herein, we cannot exhaustively investigate the variety of options and so will focus on the most obvious candidates, which are polynomials of orders 0 through 2. The zero-order polynomial is an approximation of η(x) as a uniform value over the element length (e.g., as in the Godunov conceptual model, discussed in Sect. 3) – which could be simply represented by η(x)=η_e so that the face values of η_d and η_u are ignored. Note that even this choice has alternate forms – a slightly different scheme could be constructed using $\mathit{\eta }\left(x\right)=({\mathit{\eta }}_{\mathrm{u}}+{\mathit{\eta }}_{\mathrm{d}})/\mathrm{2}$, which is also uniform but ignores η_e. A first-order polynomial for η(x) implies a linear free-surface slope across the element, which might be represented by a slope from η_d to η_u. A second-order polynomial implies a quadratic curvature to the free surface that can pass smoothly through η_u, η_e, and η_d. Clearly this idea could be extended to cubic splines by including adjacent control-volume values. Beyond these polynomials there are other options that might be suitable. For example, we could use the three discrete η values to provide piecewise linear slopes from η_u to η_e and η_e to η_d. The open-ended nature of the η(x) discretization should allow future development of a variety of finite-volume forms that can be easily demonstrated to be well-balanced and consistent with the above derivations.

Figure 7 Non-monotonic stair steps of the cross-sectional area with linear monotonic free surface. The dashed gray line is where the cross-section area reverses from increasing to decreasing.

By its definition in Eq. (37) with constraint Eq. (38), the weighting function λ(x) is an abstraction of the topographic pressure distribution over a finite volume, which is affected by both topography and the local slope of the free surface, as illustrated in Figs. 3, 5, and 6. The λ(x) is more complicated and abstract than η(x) because the free-surface elevation is approximated as uniform across the channel at any x location, but λ(x) represents the integrated effect of complex 3-D topography, e.g., the A_R(m) stair step as illustrated in Fig. 4. The key point is that the λ(x) weighting function has an integral constraint of $\int \mathit{\lambda }\mathrm{d}L=\mathrm{1}$ because the change in the cross-sectional area over the control-volume flux faces, A_d−A_u, has already been extracted from the integral of PA_R through Eq. (39). A comparison of Figs. 5 and 6 shows that λ(x) is not independent of η(x) and hence arguably should be a moderator between the static topography and the dynamics of the free surface. However, the development of λ(x, t) forms that are dynamically dependent on η(x, t) is beyond the scope of the present work. Herein, we appeal to Occam's razor and note that the simplest λ(x) that satisfies the integral constraint is $\mathit{\lambda }\left(x\right)=L_{\mathrm{e}}^{-\mathrm{1}}$, where L_e is the control-volume length, illustrated in Fig. 1. This choice of λ(x) is a zero-order polynomial implying the stair steps are uniformly distributed over the element, as illustrated in Fig. 3. Clearly, for the nonlinear cases in Figs. 5 and 6, this form of λ(x) is an approximation that reduces the nonlinear interaction between the free surface and the topography. Arguably, this is consistent with a discrete scheme using a single-valued geometry function such as η_e=f(V_e) that neglects the effect of the free-surface slope on the relationship between volume and surface elevation. Note that because the A_d−A_u area is already extracted from the integral, using $\mathit{\lambda }\left(x\right)=L_{\mathrm{e}}^{-\mathrm{1}}$ ensures that the T_e term is exact at the linear limit and a bounded approximation for nonlinear interactions. That is, with $\mathit{\lambda }\left(x\right)=L_{\mathrm{e}}^{-\mathrm{1}}$ it follows that

$${\displaystyle }{\displaystyle }{T}_{\mathrm{e}}={\displaystyle \frac{g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)}{L_{\mathrm{e}}}}\underset{L}{\int }\mathit{\eta }\left(x\right)\mathrm{d}x\le g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\max \left(\mathit{\eta }\left(x\right)\right).$$

Thus, the largest possible value for the T_e term is based on the maximum piezometric pressure over the control volume acting on the difference between upstream and downstream areas.

The simple form of $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$ is consistent with either increasing cross-sectional area (γ(x)>0 over L_e) or decreasing cross-sectional area (γ(x)<0 over L_e) but is questionable for a non-monotonic case (e.g., Fig. 6). To examine this issue in more detail, consider the somewhat simpler case of a linear free-surface slope with a cross-sectional area that increases along the streamwise direction in the upper section and decreases in the lower section, as shown in Fig. 7. The $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$ approach cannot represent the actual change in cross-section area but instead provides a linear trend between A_u and A_d as shown in Fig. 8. In contrast, if the control volume in Fig. 7 were split into two separate volumes at the centerline (i.e., the dashed gray line), then the stair steps of Fig. 7 would be readily represented by the $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$ approximation as both control volumes would be monotonic.

Figure 8 Approximated finite-element cross-section characteristics using $\mathit{\lambda }\left(x\right)=L_{\mathrm{e}}^{-\mathrm{1}}$ for non-monotonic topography in Fig. 7. The dashed gray line is where the cross-section area reverses from increasing to decreasing.

A comprehensive investigation of different forms for η(x) and λ(x) in the T_e term is clearly needed and will take significant future effort. For the present purposes, we examine the simplest polynomial forms for η(x) in combination with the zeroth-order $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$. To provide insight into how a higher-order λ(x) discretization adds complexity in the derivation, we can derive a simple linear form of λ(x) that depends only on topography and couples with a linear form of η(x). Such a λ(x) is unlikely to be useful with a dynamic free surface but is illustrative of the complexity that can be developed with quadrature of even simple linear equations. We will use the nomenclature T_e(m,n) to designate an m-order polynomial for η(x) and an n-order polynomial for λ(x).

5.2 The T_e(0,0) approximation

The simplest approximations arise by assuming uniform values such that $\mathit{\eta }\left(x\right)\approx \stackrel{\mathrm{‾}}{\mathit{\eta }}$ and $\mathit{\lambda }\left(x\right)\approx L_{\mathrm{e}}^{-\mathrm{1}}$, where $\stackrel{\mathrm{‾}}{\mathit{\eta }}$ is the average water surface elevation over the element, and we recall that ${\int }_L\mathit{\lambda }\mathrm{d}x=\mathrm{1}$. It follows that

$$\begin{array}{}\text{(50)}& {\displaystyle }{\displaystyle }{T}_{\mathrm{e}(\mathrm{0},\mathrm{0})}\equiv g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\stackrel{\mathrm{‾}}{\mathit{\eta }}.\end{array}$$

If we let ${\mathit{\eta }}_{\mathrm{e}}\approx \stackrel{\mathrm{‾}}{\mathit{\eta }}$, then Eq. (48) can be written as

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial Q_{\mathrm{e}}}{\partial t}}& {\displaystyle }-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{u}}}}-g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}+g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}\\ \text{(51)}& {\displaystyle }& {\displaystyle }=g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right){\mathit{\eta }}_{\mathrm{e}}-g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)}.\end{array}$$

Note that the RHS in the above equation is what we might infer as a discrete finite-volume version of $g\mathit{\eta }\partial A/\partial x$ in the differential source term of Eq. (47).

5.3 The T_e(1,0) approximation

If we represent the free surface by a first-order polynomial, $\mathit{\eta }\left(x\right)=ax+b$, while retaining the zeroth-order $\mathit{\lambda }\left(x\right)=L_{\mathrm{e}}^{-\mathrm{1}}$, we obtain

$$\begin{array}{}\text{(52)}& {\displaystyle }{\displaystyle }{T}_{\mathrm{e}(\mathrm{1},\mathrm{0})}=g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\left({\displaystyle \frac{aL}{\mathrm{2}}}+b\right).\end{array}$$

A linear approximation is consistent with a free surface where b=η_u and $a=-({\mathit{\eta }}_{\mathrm{u}}-{\mathit{\eta }}_{\mathrm{d}})/L$, so we obtain a finite-volume source term:

$$\begin{array}{}\text{(53)}& {\displaystyle }{\displaystyle }{T}_{\mathrm{e}(\mathrm{1},\mathrm{0})}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\left({\mathit{\eta }}_{\mathrm{u}}+{\mathit{\eta }}_{\mathrm{d}}\right).\end{array}$$

Note that the above implies products of A_dη_d and A_uη_u in the source term that can be moved to the LHS as part of the conservative flux terms. Substituting T_e(1,0) for T_e in Eq. (48) and redistributing terms provides

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial Q_{\mathrm{e}}}{\partial t}}& {\displaystyle }-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{u}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}\\ \text{(54)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g\left(A_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{u}}-A_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{d}}\right)-g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)}.\end{array}$$

Thus, the model of a linear free surface and uniform λ(x) serves to change the weighting of the gAη terms (from unity to 1∕2) and provides a source term that contains only face values of A and η, which is unlike Eq. (51) which requires an element approximation of η_e. Interestingly, the above finite-volume form does not have a differential representation as L→0. That is, the free-surface differential source term in Eq. (47) is based on $(A_{\mathrm{d}}-A_{\mathrm{u}})/L\to \mathrm{d}A/\mathrm{d}x$ as L→0. However, once we have chosen relationships for η(x) and λ(x) and moved a portion of the source term into the fluxes, an attempt to create a differential out of our source term encounters the form

$${\displaystyle }{\displaystyle }\underset{L\to \mathrm{0}}{lim}{\displaystyle \frac{A_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{u}}-A_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{d}}}{L_{\mathrm{e}}}},$$

which can only be reduced to a differential by making additional approximations as to the behavior of A and η at the faces of the finite volume. As a further insight, for the special case of a rectangular channel with frictionless flow Eq. (54) can be transformed by using A=BH, where B is the breadth and H is the depth, so that we obtain a finite-volume form of

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial }{\partial t}}\left(H_{\mathrm{e}}{U}_{\mathrm{e}}\right)& {\displaystyle }-H_{\mathrm{u}}{U}_{\mathrm{u}}^{\mathrm{2}}+H_{\mathrm{d}}{U}_{\mathrm{d}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g\left(H_{\mathrm{d}}^{\mathrm{2}}-H_{\mathrm{u}}^{\mathrm{2}}\right)\\ \text{(55)}& {\displaystyle }& {\displaystyle }=g{\displaystyle \frac{\left(H_{\mathrm{d}}+H_{\mathrm{u}}\right)}{\mathrm{2}}}\left(z_{\mathrm{u}}-z_{\mathrm{d}}\right).\end{array}$$

As L_e→0, and $S_{\mathrm{0}}=-\partial z/\partial x$, the above implies a conservative differential form of

$$\begin{array}{}\text{(56)}& {\displaystyle }{\displaystyle \frac{\partial HU}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left(H{U}^{\mathrm{2}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{\displaystyle \frac{\partial H^{\mathrm{2}}}{\partial x}}=gH{S}_{\mathrm{0}},\end{array}$$

which is commonly used in studies of simplified conservative forms (e.g., Bouchut et al., 2003; Hsu and Yeh, 2002). Thus the T_e(1,0) is also consistent with prior differential forms.

5.4 The T_e(2,0) approximation

We can take this approach further by approximating the free surface as a parabola based on {η_u, η_e, η_d}, where η_e is a characteristic free-surface height at the center of the finite volume. Using $\mathit{\eta }\left(x\right)=a{x}^{\mathrm{2}}+bx+c$ with x=0 at A_u and x=L at A_d provides

$$\begin{array}{}\text{(57)}& {\displaystyle }& {\displaystyle }c={\mathit{\eta }}_{\mathrm{u}},\text{(58)}& {\displaystyle }& {\displaystyle }b={\displaystyle \frac{\mathrm{1}}{L}}\left(-\mathrm{3}{\mathit{\eta }}_{\mathrm{u}}+\mathrm{4}{\mathit{\eta }}_{\mathrm{e}}-{\mathit{\eta }}_{\mathrm{d}}\right),\text{(59)}& {\displaystyle }& {\displaystyle }a={\displaystyle \frac{\mathrm{2}}{L^{\mathrm{2}}}}\left({\mathit{\eta }}_{\mathrm{u}}-\mathrm{2}{\mathit{\eta }}_{\mathrm{e}}+{\mathit{\eta }}_{\mathrm{d}}\right).\end{array}$$

Using $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$ results in

$$\begin{array}{}\text{(60)}& {\displaystyle }{\displaystyle }{T}_{\mathrm{e}(\mathrm{2},\mathrm{0})}={\displaystyle \frac{g}{\mathrm{6}}}\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\left({\mathit{\eta }}_{\mathrm{u}}+\mathrm{4}{\mathit{\eta }}_{\mathrm{e}}+{\mathit{\eta }}_{\mathrm{d}}\right).\end{array}$$

It is useful to multiply through and regroup terms so that A_dη_d and A_uη_u are isolated. Where these terms are balanced, they can be moved to the LHS as conservative terms. Regrouping provides

$$\begin{array}{ll}{\displaystyle }{T}_{\mathrm{e}(\mathrm{2},\mathrm{0})}& {\displaystyle }={\displaystyle \frac{g}{\mathrm{6}}}\left\{A_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}-A_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}+A_{\mathrm{d}}\left({\mathit{\eta }}_{\mathrm{u}}+\mathrm{4}{\mathit{\eta }}_{\mathrm{e}}\right)\right.\\ \text{(61)}& {\displaystyle }& {\displaystyle }\left.-A_{\mathrm{u}}\left(\mathrm{4}{\mathit{\eta }}_{\mathrm{e}}+{\mathit{\eta }}_{\mathrm{d}}\right)\right\}.\end{array}$$

So our momentum equation can be written as

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial }{\partial t}}\left(Q_{\mathrm{e}}\right)& {\displaystyle }-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{d}}}}-{\displaystyle \frac{\mathrm{5}}{\mathrm{6}}}g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}+{\displaystyle \frac{\mathrm{5}}{\mathrm{6}}}g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}g\left\{A_{\mathrm{d}}\left({\mathit{\eta }}_{\mathrm{u}}+\mathrm{4}{\mathit{\eta }}_{\mathrm{e}}\right)-A_{\mathrm{u}}\left(\mathrm{4}{\mathit{\eta }}_{\mathrm{e}}+{\mathit{\eta }}_{\mathrm{d}}\right)\right\}\\ \text{(62)}& {\displaystyle }& {\displaystyle }-g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)}.\end{array}$$

Once again, the specific representation of η(x)λ(x) provides a modification of the coefficient of the gAη terms in the conservative fluxes and sets the form of the non-conservative source term. This form does not appear to readily reduce to any differential form that is previously seen in the literature and thus provides an interesting new avenue for investigation.

5.5 The T_e(1,1) approximation

The above forms used a uniform $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$. We can readily extend the concept to analytical forms of λ(x), although it is not clear that such increasing complexity will yield an advantage in the design of a numerical model. The λ(x) function is a weighting function that reflects the distribution of the bottom elevation stair steps, as described in Eq. (4). The only restriction on λ(x) is that it must integrate to unity over L_e. Physically, as illustrated in Figs. 3, 5, and 6, the λ(x) should be a function of η(x) as well as the topography. However, we do not (as yet) have a good working framework for a dynamic representation of λ(x, η). Thus, to illustrate the complexities that arise with a nonuniform λ(x), herein we will simply analyze a somewhat arbitrary static linear relationship where λ(x)=αz(x), where z is the bottom elevation and α is a scaling constant to ensure $\underset{L}{\int }\mathit{\lambda }\mathrm{d}x=\mathrm{1}$. We introduce linear approximations of the bottom as $z\left(x\right)=Ax+B$ and the free surface as $\mathit{\eta }\left(x\right)=ax+b$. We can write these approximations as

$$\begin{array}{}\text{(63)}& {\displaystyle }& {\displaystyle }z\left(x\right)=-{\displaystyle \frac{\mathrm{1}}{L}}\left(z_{\mathrm{u}}-z_{\mathrm{d}}\right)x+z_{\mathrm{u}},\text{(64)}& {\displaystyle }& {\displaystyle }\mathit{\eta }\left(x\right)=-{\displaystyle \frac{\mathrm{1}}{L}}\left({\mathit{\eta }}_{\mathrm{u}}-{\mathit{\eta }}_{\mathrm{d}}\right)x+{\mathit{\eta }}_{\mathrm{u}}.\end{array}$$

Using $\underset{L}{\int }\mathit{\alpha }z\left(x\right)=\mathrm{1}$, it can be shown that

$$\begin{array}{}\text{(65)}& {\displaystyle }{\displaystyle }\mathit{\alpha }={\displaystyle \frac{\mathrm{2}}{L\left(z_{\mathrm{u}}+z_{\mathrm{d}}\right)}}.\end{array}$$

Using Eq. (49), a quadrature problem can be presented as

$$\begin{array}{ll}{\displaystyle \frac{T_{\mathrm{e}(\mathrm{1},\mathrm{1})}}{g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)}}=& {\displaystyle }\mathit{\alpha }\underset{L}{\int }z\left(x\right)\mathit{\eta }\left(x\right)\mathrm{d}x={\displaystyle \frac{\mathrm{2}}{L\left(z_{\mathrm{u}}+z_{\mathrm{d}}\right)}}\\ {\displaystyle }& {\displaystyle }\times \underset{L}{\int }\left[-{\displaystyle \frac{\mathrm{1}}{L}}\left(z_{\mathrm{u}}-z_{\mathrm{d}}\right)x+z_{\mathrm{u}}\right]\\ \text{(66)}& {\displaystyle }& {\displaystyle }\times \left[-{\displaystyle \frac{\mathrm{1}}{L}}\left({\mathit{\eta }}_{\mathrm{u}}-{\mathit{\eta }}_{\mathrm{d}}\right)x+{\mathit{\eta }}_{\mathrm{u}}\right]\mathrm{d}x.\end{array}$$

After some algebra, we find

$$\begin{array}{}\text{(67)}& {\displaystyle }{\displaystyle \frac{T_{\mathrm{e}(\mathrm{1},\mathrm{1})}}{g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)}}={\displaystyle \frac{\mathrm{2}{\mathit{\eta }}_{\mathrm{u}}{z}_{\mathrm{u}}+{\mathit{\eta }}_{\mathrm{u}}{z}_{\mathrm{d}}+{\mathit{\eta }}_{\mathrm{d}}{z}_{\mathrm{u}}+\mathrm{2}{\mathit{\eta }}_{\mathrm{d}}{z}_{\mathrm{d}}}{\mathrm{3}\left(z_{\mathrm{u}}+z_{\mathrm{d}}\right)}}.\end{array}$$

Unfortunately, for the above form we cannot use the redistribution trick to split T_e and move a portion to the LHS of Eq. (48). The problem is that any split term will have z_u+z_d in the denominator, which will not be conservative when used as a coefficient on a control-volume face. Thus, the T_e(1,1) form of momentum is

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial Q_{\mathrm{e}}}{\partial t}}& {\displaystyle }-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{u}}}}-g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}+g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}=g\left(A_{\mathrm{d}}-A_{\mathrm{u}}\right)\\ \text{(68)}& {\displaystyle }& {\displaystyle }\left({\displaystyle \frac{\mathrm{2}{\mathit{\eta }}_{\mathrm{u}}{z}_{\mathrm{u}}+{\mathit{\eta }}_{\mathrm{u}}{z}_{\mathrm{d}}+{\mathit{\eta }}_{\mathrm{d}}{z}_{\mathrm{u}}+\mathrm{2}{\mathit{\eta }}_{\mathrm{d}}{z}_{\mathrm{d}}}{\mathrm{3}\left(z_{\mathrm{u}}+z_{\mathrm{d}}\right)}}\right)-g{A}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)}.\end{array}$$

In general, it appears that only the uniform $\mathit{\lambda }\left(x\right)=L_{\mathrm{e}}^{-\mathrm{1}}$ form will allow us to shift part of the source term onto the flux side. Any form that has a nonuniform λ(x) must necessarily have a dependency on z(x) in the cell to satisfy $\underset{L}{\int }\mathit{\lambda }\mathrm{d}x=\mathrm{1}$, which creates terms that are inherently non-conservative. Note that this particular T_e(1,1) form is for demonstration purposes and is not recommended for use in any numerical scheme. This T_e(1,1) is predicated on an assumed weighting of λ(x)=αz(x), which does not have a physical linkage to specific cross-section geometry or expected flow conditions.

5.6 Summary of approximate forms

The T_e(0,0), T_e(1,0), T_e(2,0), and T_e(1,1) approximations all follow a similar form

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial Q_{\mathrm{e}}}{\partial t}}& {\displaystyle }-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{d}}}}-g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}\left(\mathrm{1}-{\mathit{\delta }}_{\mathrm{u}}\right)+g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}\left(\mathrm{1}-{\mathit{\delta }}_{\mathrm{d}}\right)\\ \text{(69)}& {\displaystyle }& {\displaystyle }=K_{\mathrm{e}(m,n)}-g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)},\end{array}$$

where 0≤δ_u and δ_d≤1 are fixed coefficients and K_e(m,n) is a time–space-varying topographic source term, whose exact forms are determined by the approximations used for T_e. This form was suggested in Sect. 2 with Eq. (9) based on a philosophical argument of moving as much of the RHS source terms as possible to the conservative LHS flux terms. The key point for future work is that these forms (with the exception of K_e(1,1)) are relatively straightforward in their representation of values that are the natural elements of a SVE computational model for river networks and urban drainage. This approach eliminates the need for estimating or computing the I₁ and I₂ of the Cunge–Liggett conservative form and replaces them with the simple cross-sectional area term and a K_e that is computed from discrete values of η and A. Values for δ_u, δ_d, and K_e for these forms are presented in Table 2. Examination of the above leads to the conclusion that the use of any polynomial representation of η(x) with $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$ will produce a K_e(m,0) source term that will exactly balance the piezometric pressure terms of the LHS when the free surface is flat; e.g., ${\mathit{\eta }}_{\mathrm{u}}={\mathit{\eta }}_{\mathrm{e}}={\mathit{\eta }}_{\mathrm{d}}$. Thus, these schemes are inherently well-balanced as discussed in Sect.  3. Furthermore, for these cases the source terms will be Lipschitz smooth as long as the solution variables are smooth.

Table 2 Values for conservative Aη coefficients and K_e(m,n) source terms for finite-volume schemes with (m, n) polynomial approximations of T_e. Note that forms using (m, n)=(m, 1) are not recommended without further theoretical development and are provided only for illustrative purposes.

6 Summary and discussion

The conservative differential form of the non-hydrostatic version of the Saint-Venant equations, simplified from the derivation in Sect. 4, can be written as

$$\begin{array}{ll}{\displaystyle \frac{\partial AU}{\partial t}}& {\displaystyle }+{\displaystyle \frac{\partial }{\partial x}}\left(\left[\mathit{\beta }A{U}^{\mathrm{2}}+gA\mathit{\eta }+{\displaystyle \frac{\stackrel{\mathrm{˘}}{P}}{\mathit{\rho }}}\right]\cos \mathit{\psi }\right)\\ \text{(70)}& {\displaystyle }& {\displaystyle }=g\mathit{\eta }{\displaystyle \frac{\partial A}{\partial x}}\cos \mathit{\psi }+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{A_{\mathrm{R}}}{\int }\stackrel{\mathrm{˘}}{P}\left(z\right)\mathrm{d}A-gA{S}_{\mathrm{f}}+m_{\mathrm{e}},\end{array}$$

where m_e is the source and sink of momentum from lateral fluxes per unit length (i.e., M_eL⁻¹). This equation is similar to previous work but includes both non-hydrostatic terms and effects of free-surface slope (cos ψ) that are often neglected. The key contribution of the present work is the semi-discrete, conservative, finite-volume form that corresponds to the differential form above:

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial t}}\left(U_{\mathrm{e}}{V}_{\mathrm{e}}\right)& {\displaystyle }-{\mathit{\beta }}_{\mathrm{u}}{A}_{\mathrm{u}}{U}_{\mathrm{u}}^{\mathrm{2}}\cos {\mathit{\psi }}_{\mathrm{u}}+{\mathit{\beta }}_{\mathrm{d}}{A}_{\mathrm{d}}{U}_{\mathrm{d}}^{\mathrm{2}}\cos {\mathit{\psi }}_{\mathrm{d}}-g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}\\ {\displaystyle }& {\displaystyle }+g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-{\displaystyle \frac{\cos {\mathit{\psi }}_{\mathrm{u}}}{\mathit{\rho }}}{A}_{\mathrm{u}}{\stackrel{\mathrm{˘}}{P}}_{\mathrm{u}}+{\displaystyle \frac{\cos {\mathit{\psi }}_{\mathrm{d}}}{\mathit{\rho }}}{A}_{\mathrm{d}}{\stackrel{\mathrm{˘}}{P}}_{\mathrm{d}}\\ {\displaystyle }& {\displaystyle }=g\left(A_{\mathrm{d}}\cos {\mathit{\psi }}_{\mathrm{d}}-A_{\mathrm{u}}\cos {\mathit{\psi }}_{\mathrm{u}}\right)\underset{L}{\int }\mathit{\eta }\left(x\right)\mathit{\lambda }\left(x\right)\mathrm{d}x\\ \text{(71)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\underset{L}{\int }\mathit{\gamma }\left(x\right)\underset{A_{\mathrm{R}}}{\int }\stackrel{\mathrm{˘}}{P}(x,y,z)\mathrm{d}A\mathrm{d}x-g{V}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)}+M_{\mathrm{e}},\end{array}$$

where ${\stackrel{\mathrm{˘}}{P}}_{\mathrm{u}}$ and ${\stackrel{\mathrm{˘}}{P}}_{\mathrm{d}}$ are the average non-hydrostatic pressures on the upstream and downstream cross-sectional areas. In this finite-volume form, the only approximations introduced are listed as follows: (1) uniform-density incompressibility, (2) the effects of momentum redirection around bends is either negligible or is handled in friction terms, (3) the cross-channel variability in the free-surface slope is negligible, and (4) a friction slope model can be used to represent integrated viscous effects. In addition we have introduced a geometric restriction that the upstream and downstream cross sections must be vertical planes that are orthogonal (in the horizontal plane) to the mean flow direction. The above form can be used to analyze systems that include non-hydrostatic pressure and slope gradients beyond the small-slope approximation of the traditional SVEs. As warning for future development of non-hydrostatic methods, note that the fundamental 1-D derivation is effectively treating the non-hydrostatic pressure gradients in the horizontal as absorbing and redirecting momentum along the curving channel. Thus a part of this term is, in effect, encapsulated within the approximation that allows A_u and A_d to be cross sections that are not parallel.

The principal feature of the new finite-volume formulation is the topographic source term ∫η(x)λ(x)dx that can be represented by analytical functions to approximate a smoothly varying free surface and its interaction with topography across the finite-volume element. Discrete polynomial representations of ∫η(x)λ(x)dx are evaluated in Sect. 5, with the resulting topographic terms designated as T_e(m,n) for an m-degree polynomial representing η and an n-degree polynomial representing λ. In an approximate form, the T_e(m,n) term is split into a δ_u(m,n) factor applied to gA_uη_u and a δ_d(m,n) factor applied to gA_dη_d, which become part of the conservative flux terms. The remainder of the T_e(m,n) becomes a K_e(m,n) source term in the approximate finite-volume form. Simpler conservative finite-volume forms that use the common approximations of the hydrostatic equations ($\stackrel{\mathrm{˘}}{P}\approx \mathrm{0}$) with small free-surface slope (cos ψ≈1), uniform cross-sectional velocity (β≈1), and no momentum sources (M_e=0) can be written as

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial Q_{\mathrm{e}}}{\partial t}}& {\displaystyle }-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{d}}}}-g{A}_{\mathrm{u}}{\mathit{\eta }}_{\mathrm{u}}\left(\mathrm{1}-{\mathit{\delta }}_{\mathrm{u}(m,n)}\right)\\ \text{(72)}& {\displaystyle }& {\displaystyle }+g{A}_{\mathrm{d}}{\mathit{\eta }}_{\mathrm{d}}\left(\mathrm{1}-{\mathit{\delta }}_{\mathrm{d}(m,n)}\right)=K_{\mathrm{e}(m,n)}-g{L}_{\mathrm{e}}{A}_{\mathrm{e}}{S}_{\mathrm{f}\left(\mathrm{e}\right)},\end{array}$$

where the discrete K_e and δ terms are shown in Table 2. It is worthwhile to compare the above to a finite-volume form derived using the Cunge–Liggett form of the SVEs, which could be written (using the I₁ and I₂ definitions) as

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{e}}{\displaystyle \frac{\partial Q}{\partial t}}& {\displaystyle }-{\displaystyle \frac{Q_{\mathrm{u}}^{\mathrm{2}}}{A_{\mathrm{u}}}}+{\displaystyle \frac{Q_{\mathrm{d}}^{\mathrm{2}}}{A_{\mathrm{d}}}}-g\underset{H\left(u\right)}{\int }(H-z)B\left(z\right)\mathrm{d}z\\ {\displaystyle }& {\displaystyle }+g\underset{H\left(d\right)}{\int }(H-z)B\left(z\right)\mathrm{d}z=g\underset{L}{\int }\underset{H}{\int }(H-z){\displaystyle \frac{\partial B}{\partial x}}\mathrm{d}x\\ \text{(73)}& {\displaystyle }& {\displaystyle }+g\underset{L}{\int }A\left(S_{\mathrm{0}}-S_{\mathrm{f}}\right)\mathrm{d}x.\end{array}$$

Performance of the traditional scheme depends on the specification for B(x, z) that defines the irregular bathymetry of the channel. Although the B(x, z) term is developed without any approximations, it is a nontrivial matter to simplify these terms to create practical computational forms for irregular cross-section geometry. The source term on the RHS of the Cunge–Liggett form is effectively an integration of both variations in the channel topography and the water surface elevation over the volume – similar to the new T_e – but without a limiting constraint (i.e., $\partial B/\partial x$ is not inherently limited in magnitude for irregular topography). Furthermore, the selection of B(x, z) for the RHS affects the integration on the LHS hydrostatic pressure term – thus obtaining a well-balanced source term that compensates for S₀ will also affect the conservative flux terms. The λ(x) used to develop the K_e term Eq. (72) serves a similar purpose to the source-term integral in the Cunge–Liggett form, but it provides a simple weighting function that can be analytically integrated with an approximation for η(x) and is inherently constrained such that $\int \mathit{\lambda }\mathrm{d}x=\mathrm{1}$ over a control volume. The other major difference between the two approaches is that Eq. (73) uses S₀ as a source term on the RHS of the equations, whereas in the new approach Eq. (72) dispenses with this artifice so that the source terms only include friction (which is guaranteed to damp momentum) and the portion of the topographic effects that cannot be transferred into the conservative δ_u and δ_d terms.

There is a long history of the use of S₀ in the source term in the SVEs, and it is indeed part of the author's prior model (Liu and Hodges, 2014). However, the use of S₀ with irregular geometry brings the problems of creating a well-balanced conservative scheme, as discussed in Sect. 3. Furthermore, the use of S₀ in the source term requires the pressure term to be treated as the hydrostatic pressure rather than the piezometric pressure, as shown in Sect. 4. Because the hydrostatic pressure is a function of depth, its integration over a cross-sectional area requires knowledge of the distribution of depth across the channel – a significant computational complexity for irregular topography. In contrast, the integration of the piezometric pressure over a cross section is exactly PA and does not require knowledge of how depth is distributed across the channel. Other authors have noted similar problems: Schippa and Pavan (2008) derived a conservative differential form that retained gI₁ in the flux terms and removed S₀ by showing that it could be combined with the I₂ source term as $\partial I_{\mathrm{1}}/\partial x$ for a uniform water level. It seems that the Schippa and Pavan (2008) differential equation might be preferred to the approach proposed herein for high-resolution topography with small grid spacing (i.e., where we have confidence that the computation of I₁ is meaningful). However, at larger scales where geometric cross sections are broadly spaced and the computation of I₁ is questionable, the simplicity of using A and η for piezometric pressure gradient terms is likely to be preferred. Other authors, notably Rosatti et al. (2011), have simply accepted $gA\partial \mathit{\eta }/\partial x$ as an unavoidable source term rather than dealing with the problems of obtaining a well-balanced method with the hydrostatic pressure.

Since S₀ was not in de Saint-Venant's original paper, how did it come to be commonly used in the Saint-Venant equations? Arguably there are two sources associated with different simulation scales: (1) in hydrology the kinematic wave model provides S₀=S_f, which leads to prioritization of S₀ as a hydraulic parameter; and (2) in mathematics the equation

$$\begin{array}{}\text{(74)}& {\displaystyle }{\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial hu}{\partial x}}=s\end{array}$$

is the canonical form of a 1-D inhomogeneous hyperbolic advection equation for depth h and velocity u, which leads to a prioritization of the depth as a fundamental parameter and requires that S₀(x) be relegated to the source term for irregular geometry. However, for solution of the SVEs there is no need to exactly mimic the hyperbolic advection equation to obtain a conservative form; thus there is no need to introduce S₀. The above comparisons of the Cunge–Liggett form and the new finite-volume form illustrate the additional complexity of introducing S₀. Beyond these issues is a more fundamental problem: numerical methods for inhomogeneous partial differential equations are only well-posed if the source terms are Lipschitz smooth (e.g., Iserles, 1996) – otherwise one should not be surprised by numerical instabilities and/or difficulties in convergence. Any river model that uses raw data from surveyed cross sections will inherently have a non-smooth S₀(x). As a result, much of the computational complexity is likely to be compensating for the lack of Lipschitz smoothness in the boundary conditions of topography. In contrast, when the free-surface elevation (piezometric pressure) is used instead of depth (hydrostatic pressure), then S₀(x) disappears from the SVEs and the smoothness of the source term is assured (for smooth solution variables) by the choice of smooth functions for λ(x) and η(x) and the friction slope model S_f(x). In general, as long as the solutions of η and Q remain smooth the source term of the SVEs, as derived herein, should remain smooth. That is, the approach herein cannot guarantee a smooth solution, but it can guarantee that any observed non-smoothness in the source terms during a simulation is a result of non-smoothness in the solution variables rather than direct forcing of boundary conditions.

As a matter of pure speculation, the new form of the finite-volume equations brings up some interesting possibilities for large-scale modeling. Imagine that we would like to model a river network or urban drainage network where we have some high-resolution (1×1 m) data in some areas but not in others. Let us also say our computer power limits us to a SVE solution with a median cell length of about 20 m. Can we use our high-resolution knowledge directly at the coarse scale? Arguably, the λ(x) approach could give us a means of directly incorporating effects of subgrid-scale topography into a single source term – however, significant theoretical work is required to develop a consistent methodology that retains the fundamental well-balanced and Lipschitz smooth characteristics of the finite-volume equations.

Clearly, there remains much practical experimentation to be done in comparing various forms of T_e(m,n) and the effectiveness of different numerical solution methods with different η(x) polynomial representations. There is also a need to develop a firm theoretical relationship between λ(x) and η(x) to overcome the difficulties illustrated with Figs. 5 and 6 where nonlinear interactions between the free surface and topography cannot be readily represented with the uniform $\mathit{\lambda }\left(x\right)=L^{-\mathrm{1}}$ applied herein. Finally, there are substantial questions on whether the new approach can be used with methods designed to satisfy the entropy criterion (e.g., Greenberg and Leroux, 1996; Harten et al., 1983; Lax, 1973). It is hoped that researchers will consider adapting their finite-volume codes to the form of Eq. (72) and reporting their experience.

7 Conclusions

New finite-volume forms of integral momentum equations for unsteady flow in open channels with varying cross sections are derived in this paper. These equations reduce to the classic differential forms of the Saint-Venant equations (under the T_e(0,0) and T_e(1,0) approximations) and also provide new approximate finite-volume forms that are suitable for analytical representations of topography and free-surface elevation over a finite-volume element. The new forms use the piezometric pressure (free-surface elevation) rather than hydrostatic pressure (depth) as a fundamental variable and thus do not include the channel slope, S₀. This approach provides a cleaner finite-volume form as the nonlinear interactions of topography and the free surface are handled in a single integral term, $g(A_{\mathrm{d}}-A_{\mathrm{u}})\int \mathit{\eta }\mathit{\lambda }\mathrm{d}x$, where λ is a quadrature weighting function and A_d−A_u is the downstream increase in cross-sectional area of the control volume. The introduction of λ(x) provides a potential avenue to convert practical knowledge of fluid and/or geometry variations within a control volume into source and flux terms of the finite-volume equations. The derivations herein can be used to generate a variety of conservative, well-balanced, finite-volume forms for the SVEs that could be employed with a wide range of numerical discretization schemes. This work provides only the theoretical development for the new finite-volume equations; the practical implementation and numerical testing remains a subject for future work.

Appendix A

To elaborate on Eq. (19), the advection of some quantity ϕ through a control volume bounded by surface S can be represented as

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }\underset{S}{\oint }\mathit{\varphi }{u}_k{n}_k\mathrm{d}A=\underset{S}{\oint }\mathit{\varphi }{u}_{\widehat{\mathit{n}}}\mathrm{d}A,\end{array}$$

where $u_{\widehat{\mathit{n}}}$ is the local normal velocity across a flux surface. For a flux surface of area A, we define the nonlinearity of distribution across the surface as

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\displaystyle }\mathit{\alpha }\equiv {\displaystyle \frac{\mathrm{1}}{A\mathrm{\Phi }{U}_{\widehat{\mathit{n}}}}}\underset{A}{\int }\mathit{\varphi }{u}_{\widehat{\mathit{n}}}\mathrm{d}A,\end{array}$$

where Φ is the average of ϕ over flux area A, and $U_{\widehat{\mathit{n}}}$ is the average surface normal velocity. Note that because $u_{\widehat{\mathit{n}}}$ and $U_{\widehat{\mathit{n}}}$ are derived from projection of the velocity onto the surface normal vector, they have positive signs for outflows and negative sign for inflows. Consistency of $u_{\widehat{\mathit{n}}}$ and $U_{\widehat{\mathit{n}}}$ ensures that α≥0; however, this introduces a nomenclature difficulty because the sign does not depend solely on the nominal upstream or downstream direction of the flow. To simplify the nomenclature, we consider only upstream (u) and downstream (d) flux surfaces for a control volume and represent the normal velocity by U_⟂ such that a downstream velocity on any face is $U_{⟂}>\mathrm{0}$ and an upstream velocity on any face is $U_{⟂}<\mathrm{0}$. It follows that

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\displaystyle }\underset{S}{\oint }\mathit{\varphi }{u}_k{n}_k\mathrm{d}A=-{\left(\mathit{\alpha }\mathrm{\Phi }{U}_{⟂}A\right)}_{\mathrm{u}}+{\left(\mathit{\alpha }\mathrm{\Phi }{U}_{⟂}A\right)}_{\mathrm{d}},\end{array}$$

where

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\displaystyle }\mathit{\alpha }\equiv {\displaystyle \frac{\mathrm{1}}{A\mathrm{\Phi }{U}_{⟂}}}\underset{A}{\int }\mathit{\varphi }{u}_{⟂}\mathrm{d}A.\end{array}$$

Let $Q\equiv U_{⟂}A$ with the same sign conventions as U_⟂, so that

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\displaystyle }\underset{S}{\oint }\mathit{\varphi }{u}_k{n}_k\mathrm{d}A=-(\mathit{\alpha }\mathrm{\Phi }Q{)}_{\mathrm{u}}+(\mathit{\alpha }\mathrm{\Phi }Q{)}_{\mathrm{d}}.\end{array}$$

Consider a channel with a linear free surface so that a velocity parallel to the free surface ($u_{\widehat{\mathit{i}}}$) has the same direction over the entire control volume. Let $\mathit{\varphi }=u_{\widehat{\mathit{i}}}$ and $\mathrm{\Phi }=U_{\widehat{\mathit{i}}}$ so we obtain

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\displaystyle }\underset{S}{\oint }{u}_{\widehat{\mathit{i}}}{u}_k{n}_k\mathrm{d}A=-{\left(\mathit{\alpha }{U}_{\widehat{\mathit{i}}}Q\right)}_{\mathrm{u}}+{\left(\mathit{\alpha }{U}_{\widehat{\mathit{i}}}Q\right)}_{\mathrm{d}},\end{array}$$

where $U_{\widehat{\mathit{i}}}$ is interpreted as the average velocity parallel to the free surface with the same sign conventions as U_⟂. For a free surface at angle ψ from the horizontal it follows that

$$\begin{array}{}\text{(A7)}& {\displaystyle }{\displaystyle }\cos \mathit{\psi }={\displaystyle \frac{U_{⟂}}{U_{\widehat{\mathit{i}}}}}.\end{array}$$

Using Eq. (A4) with $\mathrm{\Phi }=U_{\widehat{\mathit{i}}}$ and $U_{⟂}=U_{\widehat{\mathit{i}}}\cos \mathit{\psi }$ provides

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\displaystyle }\mathit{\alpha }={\displaystyle \frac{\mathrm{1}}{A{\left[U_{\widehat{\mathit{i}}}\right]}^{\mathrm{2}}\cos \mathit{\psi }}}\underset{A}{\int }{u}_{\widehat{\mathit{i}}}{u}_{⟂}\mathrm{d}A.\end{array}$$

Noting that $u_{\widehat{\mathit{i}}}=u_{⟂}{\cos }^{-\mathrm{1}}\mathit{\psi }$ we obtain

$$\begin{array}{}\text{(A9)}& {\displaystyle }{\displaystyle }\mathit{\alpha }={\displaystyle \frac{\mathrm{1}}{A{\left[U_{\widehat{\mathit{i}}}\right]}^{\mathrm{2}}}}\int {\left[u_{\widehat{\mathit{i}}}\right]}^{\mathrm{2}}\mathrm{d}A.\end{array}$$

Thus, α of Eq. (A4) for $\mathrm{\Phi }=U_{\widehat{\mathit{i}}}$ is identical to β of Eq. (18). It follows that

$$\begin{array}{}\text{(A10)}& {\displaystyle }{\displaystyle }\underset{S}{\oint }{u}_{\widehat{\mathit{i}}}{u}_k{n}_k\mathrm{d}A=-{\left(\mathit{\beta }{U}_{\widehat{\mathit{i}}}Q\right)}_{\mathrm{u}}+{\left(\mathit{\beta }{U}_{\widehat{\mathit{i}}}Q\right)}_{\mathrm{d}},\end{array}$$

which is Eq. (19). However, because $Q\equiv U_{⟂}A$ it follows that $Q=U_{\widehat{\mathit{i}}}A\cos \mathit{\psi }$. Thus, equivalent forms are

$$\begin{array}{}\text{(A11)}& {\displaystyle }\underset{S}{\oint }{u}_{\widehat{\mathit{i}}}{u}_k{n}_k\mathrm{d}A& {\displaystyle }=-{\left(\mathit{\beta }{U}_{\widehat{\mathit{i}}}^{\mathrm{2}}A\cos \mathit{\psi }\right)}_{\mathrm{u}}+{\left(\mathit{\beta }{U}_{\widehat{\mathit{i}}}^{\mathrm{2}}A\cos \mathit{\psi }\right)}_{\mathrm{d}},\text{(A12)}& {\displaystyle }& {\displaystyle }=-{\left(\mathit{\beta }{\displaystyle \frac{Q^{\mathrm{2}}}{A\cos \mathit{\psi }}}\right)}_{\mathrm{u}}+{\left(\mathit{\beta }{\displaystyle \frac{Q^{\mathrm{2}}}{A\cos \mathit{\psi }}}\right)}_{\mathrm{d}},\end{array}$$

which all reduce to conventional forms with Q=UA when cos ψ=1.