The Saint-Venant equations are commonly used as the governing equations to solve for modeling the spatially varied unsteady flow in open channels. The presence of uncertainties in the channel or flow parameters renders these equations stochastic, thus requiring their solution in a stochastic framework in order to quantify the ensemble behavior and the variability of the process. While the Monte Carlo approach can be used for such a solution, its computational expense and its large number of simulations act to its disadvantage. This study proposes, explains, and derives a new methodology for solving the stochastic Saint-Venant equations in only one shot, without the need for a large number of simulations. The proposed methodology is derived by developing the nonlocal Lagrangian–Eulerian Fokker–Planck equation of the characteristic form of the stochastic Saint-Venant equations for an open-channel flow process, with an uncertain roughness coefficient. A numerical method for its solution is subsequently devised. The application and validation of this methodology are provided in a companion paper, in which the statistical results computed by the proposed methodology are compared against the results obtained by the Monte Carlo approach.

Unsteady open-channel flows are a common occurrence in hydrology and hydraulics problems. They arise as a result of the movement of water waves in natural or artificial channels (Sturm, 2001). Understanding and tracing the movement of such water waves along the channels is of great importance in addressing engineering flow problems, including flood forecasting, flood control, hydrograph generation, and several others (Chow, 1959). The technique used to approximate and trace such water waves is known as flood routing, and the governing equations that are commonly used to solve for the unsteady flows in flood routing problems are known as the Saint-Venant equations (Chanson, 2004).

Various uncertainties may add to the complexity of solving the Saint-Venant equations (Gates and AlZahrani, 1996a; Ercan and Kavvas, 2012a) and these may correspond to several factors. Physical conditions of open channels may be uncertain due to their high degree of variability (Sturm, 2001). One example is Manning's roughness coefficient, which greatly depends on the channel vegetation, bed material, bedforms, and even on the position of the free water surface (Chow, 1959; Sturm, 2001; Ercan and Kavvas, 2012a). With the uncertainties in quantifying or characterizing these factors, the roughness coefficient becomes extremely difficult to estimate (Sturm, 2001), rendering it uncertain. Channel geometric parameters may also be uncertain. This includes the channel bed slope (Ercan and Kavvas, 2012a) and the channel cross section geometry, the latter of which may exhibit significant spatial variability across a river due to its irregular form and due to the changes it may undergo along the direction of flow (Chow, 1959). Other uncertainties may also arise from lateral inflows and initial conditions due to their spatial and/or temporal variability (Liang and Kavvas, 2008), as well as from the upstream boundary conditions due to the temporal variability of the inflows into the channel.

As a result of such uncertainties, the channel and flow parameters may be considered spatially and/or temporally random at the local scale of a river cross section, rendering the system behavior uncertain (Gates and AlZahrani, 1996a). Therefore, deterministically solving the Saint-Venant equations in this case would no longer be providing a representative solution to the flood routing problem being considered. In fact, in this study the governing partial differential equations (PDEs), i.e., the Saint-Venant equations, will be transformed into stochastic PDEs because the channel and flow parameters are now stochastic and can be described as random functions (Liang and Kavvas, 2008). This means that the dependent variables that will be solved for by these equations (e.g., flow velocity and depth) will also be spatiotemporal random functions. Hence, instead of solving for the deterministic values of the dependent variables, the goal will be to solve for their statistical properties (Van Kampen, 1976), which can be obtained at designated discrete time–space positions.

Two popular methods can be used for such a stochastic solution to nonlinear problems: the finite-order analysis and the Monte Carlo (MC) approach. Applying the expectation operator, or any other statistical moment operator, to nonlinear difference equations with stochastic parameters may result in nonlinear expressions that are difficult, or even impossible, to simplify into terms involving the known moments of the model parameters and the unknown moments of the dependent variables (Gates and AlZahrani, 1996a). Finite-order analysis overcomes this by performing a Taylor series expansion of the difference equations about the expected values of the parameters, from which higher-order terms are truncated. For example, truncating the Taylor series after the first-order term is known as the first-order approximation, which is a good approximation when the system nonlinearity is not too high, and when the stochastic parameters have relatively small coefficients of variation (Dettinger and Wilson, 1981). However, with highly nonlinear problems, instead of using higher-order approximations of the finite-order method, it may be required and more efficient to use full-distribution methods such as the MC approach (Dettinger and Wilson, 1981).

The MC approach is well-known for simulating differential equations with stochastic parameters, and is used to determine the distributions of the unknown stochastic dependent variables (Freeze, 1975; Smith and Freeze, 1979; Bellin et al., 1992). This method involves repeatedly solving the governing equations in a deterministic fashion, varying the stochastic parameters for each run, in order to obtain a set of several realizations for each of the dependent variables. When a sufficient number of realizations is obtained, they can be used to determine the required statistical properties, including the mean system behavior and the standard deviation (Gates and AlZahrani, 1996b). Therefore, the MC simulations require two models: one which generates realizations for the stochastic parameters, and another (finite-difference model) which deterministically solves the governing flow equations for each realization (Gates and AlZahrani, 1996b). The MC approach is generally accepted as the most robust approach for uncertainty evaluation, as well as the benchmark for comparing other new methods (Scharffenberg and Kavvas, 2011). The full distribution characteristics may be estimated using the MC approach, which is more intuitive than the finite-order methods (Gates and AlZahrani, 1996a). However, the main drawback of the MC approach is its computational expense due to the usual running of a large number of simulations of the process under study in order to obtain accurate results (Dettinger and Wilson, 1981).

To bypass the need for solving the unsteady open-channel flow governing equations several times, a new methodology is proposed in this study in order to solve for the expected system behavior and variability in only one simulation. This methodology involves upscaling the governing stochastic differential equations from the point scale (at which they are originally valid) to the field scale. Ensemble averaging has been a common approach to upscale hydrologic equations that are linear (Gelhar and Axness, 1983; Kitanidis, 1988; Rubin and Dagan, 1989; Kapoor and Gelhar, 1994; Kavvas and Karakas, 1996; Wood and Kavvas, 1999a, b) or nonlinear (Mantoglou and Gelhar, 1987; Tayfur and Kavvas, 1994; Horne and Kavvas, 1997; Dogrul et al., 1998), in which case these equations are averaged to become deterministic differential equations. These developed deterministic differential equations use statistical descriptions, such as the mean and variance, to represent the values of the stochastic parameters (Liang and Kavvas, 2008). However, most of the studies performing the ensemble averaging technique on nonlinear conservation equations used the regular perturbation method, which includes linearization assumptions and which only works for small fluctuations in the dependent variables (Kavvas, 2003). Other techniques that have been applied, which are not limited by small fluctuations, include the decomposition method (Serrano, 1995), a combination of volume averaging with nonlinear dynamics (Duffy, 1996; Duffy and Cusumano, 1998), as well as the theory of fractals and multifractals (Puente, 1996). Nonetheless, due to some limitations of such methods when used for stochastic nonlinear hydrologic processes, the upscaling method used in this study is chosen to be that of Kavvas (2003).

Kavvas (2003) developed general ensemble average conservation equations (to second order) for nonlinear and linear hydrologic processes in order to determine their probabilistic and mean behavior. The “master key” equations developed may be used on any stochastic hydrologic process after being rewritten as one or more linear/nonlinear stochastic ordinary differential equations (ODEs). This utilization leads to a special Lagrangian–Eulerian form of the Fokker–Planck equation (LEFPE) that models the time–space evolution of the probability density of the dependent variables of any nonlinear/linear stochastic dynamic process (Kavvas, 2003). Such a methodology has been successfully applied to many hydrologic processes, including unsaturated water flow (Kim et al., 2005b), root-water uptake (Kim et al., 2005a), solute transport (Liang and Kavvas, 2008), snow accumulation and melt (Ohara et al., 2008), unconfined groundwater flow (Cayar and Kavvas, 2009a, b), and kinematic open-channel flow (Ercan and Kavvas, 2012a, b).

Noting that the characteristic forms of the Saint-Venant equations are nonlinear ODEs, it is proposed in this study to apply to them their corresponding master key equation from Kavvas (2003). From this operation the corresponding LEFPE of the Saint-Venant open-channel flow equations is obtained, thus providing the ability to model the uncertainties of the channel and flow parameters and to compute their effect on the behavior of the system. Therefore, under the appropriate initial and boundary conditions, the probability density functions (PDFs) of the dependent variables can be computed (to exact second order) through the LEFPE, and the ensemble behavior of the system can be described.

The advantages of using the LEFPE in tackling the flood routing problem greatly echo those of the classical Fokker–Planck equation (FPE). In fact, the LEFPE directly solves for the PDFs of the dependent variables of the system in both time and space, it is linear in the variable being solved for (i.e., in the PDF), and unlike the many simulations usually performed for the MC approach, the LEFPE produces the complete ensemble model results with only one single simulation. As such, the LEFPE provides not only the mean and variance of the process but also a complete description of the evolution of the dependent variables' PDFs in a computationally efficient manner. Note that the LEFPE does not make any linearization assumptions, it works with a wide-ranging parameter space, and the only assumption about the physical process it makes is the finite correlation time for the process (Ercan and Kavvas, 2012a).

Therefore, following from the above discussion, the main objective of this study is to apply the upscaling method based on the LEFPE approach in Kavvas (2003) to the characteristic form of the stochastic Saint-Venant equations in order to derive a new methodology that solves for the probability density of the dependent flow variables, and that quantifies the expected behavior and variability of the system in one shot, instead of running a large number of simulations.

The Saint-Venant equations, also known as the spatially varied unsteady flow equations (Sturm, 2001), are the two governing equations used to describe an unsteady open-channel flow problem that will be solved using the hydraulic routing technique (Chow, 1959; Viessman et al., 1977; Sturm, 2001). They consist of the continuity equation and the momentum equation which are used simultaneously in order to solve for the two unknowns (velocity and depth, or discharge and depth). The naming of these equations comes from the French mathematician Adhémar-Jean-Claude Barré de Saint-Venant who published the equations describing one-dimensional unsteady open-channel flow in 1871 (Barré de Saint-Venant, 1871).

Several assumptions are made when deriving these equations
(Viessman et al., 1977; Sturm, 2001), including
unidirectional flow and uniform cross-sectional velocity, hydrostatic
pressure, small channel bed slope, steady state estimation of friction loss,
and incompressible flow. Following these assumptions, the Saint-Venant
equations for unsteady open-channel flow of an incompressible fluid in a
rectangular, prismatic channel (with no lateral inflow/outflow) may be
written as follows (Viessman et al., 1977):

Since closed-form solutions to the Saint-Venant equations have not been
obtained due to the presence of nonlinear terms, it has not been possible to
solve these equations analytically except when extreme simplifications are
applied (Sturm, 2001; Chaudhry, 2008). As a result,
several numerical techniques have been developed in order to solve the
Saint-Venant equations deterministically in their full form, without major
simplifications. The most frequently used of these techniques are
finite-difference methods (Abbott and Ionescu, 1967; Fread, 1973; Beam
and Warming, 1976; Fennema and Chaudhry, 1986; Garcia and Kahawita, 1986;
Venutelli, 2002), which solve the governing equations explicitly or
implicitly along a fixed or adaptive

Furthermore, another approach to solve the Saint-Venant equations may be
utilized after realizing that these equations are hyperbolic PDEs
(Chaudhry, 2008). This approach is known as the method of
characteristics (MOC; Abbott, 1966), which is one of the earliest
and most exact methods for solving hyperbolic PDEs (Tannehill
et al., 1997), and which was used early on for solving the Saint-Venant
equations (Amein, 1966; Woolhiser and Liggett, 1967; Lai, 1988). The MOC
may be used to transform a hyperbolic PDE into a system of ODEs, which may
be simpler to solve (Sturm, 2001). These ODEs are usually
divided into two equations: the characteristic equation (i.e., the ODE
describing the characteristic path), and the compatibility equation (i.e.,
the ODE that describes the process behavior along that characteristic path)
(Hoffman, 2001). After this transformation by the MOC, the
finite-difference approximations of the derivatives can then be applied to
the characteristic form of the hyperbolic PDE, instead of applying them to
its original form. Note that in time and one-space dimensions, the
characteristic equations represent curves in the

From the several techniques available to solve for the Saint-Venant
equations, the MOC is chosen for this study. This is because, as was
mentioned in Sect. 1, the upscaling technique based on the LEFPE approach in
Kavvas (2003) can be applied to hydrologic processes that are
written as one or more ODEs. As such, it is imperative for the progression
of this study to transform the Saint-Venant equations into their
characteristic form in order to write them as a system of ODEs. With two
characteristic directions, the Saint-Venant equations are transformed by the
MOC into a system of four ODEs: two characteristic equations and two
corresponding compatibility equations. When finite-difference approximations
are applied to the characteristic form of the Saint-Venant equations, the
results can be numerically computed along an irregular

Through a linear combination of the continuity and momentum equations (Eqs. 1
and 2), the characteristic equations for unsteady open-channel flow of an incompressible
fluid in a rectangular, prismatic channel with no lateral inflow can be
written as follows (Sturm, 2001):

Equations (3) and (5) represent two different velocity expressions defining the two characteristic
directions of the Saint-Venant equations: the former defining the positive
characteristic curve (

Therefore, the MOC transforms the two governing PDEs into a system of four ODEs that are differentiated with respect to time only, and no longer with respect to space. This transformation provides the ability to use the upscaling technique based on the LEFPE approach in Kavvas (2003) on the Saint-Venant equations in order to derive this study's proposed methodology involving the ensemble-averaged equations of stochastic unsteady open-channel flow.

In this section, a new methodology for solving the stochastic Saint-Venant equations will be introduced and derived, and a numerical discretization scheme will be devised for it as well. The proposed methodology aims at obtaining the statistical properties of the dependent variables of the unsteady open-channel flow system in only one simulation, as opposed to the large number of simulations usually involved in the MC approach. The following derivation involves assuming the Manning's roughness coefficient as an uncertain parameter, but similar steps can be followed even when the uncertainty is assumed to arise from other parameters.

Following Kavvas (2003), a system of point-scale conservation equations can be
written for a dynamical system as follows:

Solving the LEFPE, Eq. (9), under the appropriate
initial and boundary conditions provides the spatiotemporal evolution of the
PDF of the vector of state variables (

Since the LEFPE was developed for a system of ODEs (Eq. 7), and since the characteristic form of the
Saint-Venant equations is a system of four nonlinear ODEs (Eqs. 3 to 6), it is
proposed to apply to these equations the corresponding LEFPE after making
some substitutions and adjustments. First, the friction slope (

Note that in equations to follow, the

Note that the LEFPE has the form of an advection–diffusion equation. In
Eq. (19), the first four terms on the right-hand side
represent the advection terms, while the remaining terms represent the
diffusion terms. Within the advection terms, the expected values of the

In order to apply the derived FPE methodology, the FPE represented in Eq. (26) must be solved using an appropriate numerical scheme. In general, finite-difference and finite-element methods, among others, have been widely used to solve FPEs numerically. Many studies have compared several such methods to determine how they perform against each other in solving different FPEs. In their study, Pichler et al. (2011) looked at the central finite-difference method, the alternating directions implicit (ADI) method, as well as finite-element methods. They mentioned that finite-difference methods are computationally more economical than finite-element methods as the number of dimensions increases. Park and Petrosian (1996) performed a comparison between the Chang and Cooper (1970) scheme, the schemes presented in Larsen et al. (1985), and some implicit schemes (including a fully implicit mid-point difference method). They also studied the semi-implicit forms of these schemes (e.g., the Crank–Nicholson scheme). They concluded that, among these schemes, the best finite-difference method for solving their FPEs was the Chang–Cooper scheme, as it was the most robust and most stable over a wide range of parameters among the other methods tested in their study.

In fact, the Chang–Cooper scheme has been cited as one of the most widely known schemes for solving the classical FPE numerically. In their paper, Chang and Cooper (1970) developed a practical numerical differencing scheme for the solution of the one-dimensional classical FPE. This scheme uses a centered difference for the diffusion term, a weighted difference for the advection term, and requires that the quasi-equilibrium solution to the FPE be satisfied exactly at any given time along the mesh nodes. The Chang-Cooper scheme, which has first-order convergence in space and time, is a widely used scheme that ensures the non-negativity of the solution, the conservation of the probability mass (in the absence of any external sources or sinks), and the exact representation of the analytical solution upon equilibration. As a result, the Chang–Cooper scheme is highly accurate with a relatively small number of required mesh nodes. While Chang and Cooper (1970) developed and applied that scheme to a one-dimensional FPE, Kim et al. (2005a) generalized and applied the Chang–Cooper scheme to a two-dimensional FPE. In a similar manner, for the case of this study, an attempt is made to generalize the Chang–Cooper scheme to the four-dimensional FPE shown in Eq. (26).

First, Eq. (26) is rewritten as follows:

This study proposed a new methodology to model the expected behavior and variability of a system described by the stochastic open-channel flow equations. The governing equations that were used to represent the flood routing problem in this study are the continuity and momentum equations, otherwise known as the Saint-Venant equations. Many uncertainties can add to the complexity of solving the Saint-Venant equations in engineering routing problems. These uncertainties may include uncertainties in the channel's physical and geometric properties, as well as uncertainties in the lateral inflows and upstream boundary conditions, all of which render the Saint-Venant equations stochastic. As such, the dependent variables that will be solved for by these equations will also become stochastic, thus requiring that their statistical properties be solved for at specific time–space locations. Therefore, with uncertain parameters, the Saint-Venant equations have to be solved within a stochastic framework in order to quantify the ensemble behavior and variability of the system being considered. While the Mote Carlo method is a viable approach for the solution of such a stochastic unsteady open-channel flow problem, its computational expense and its large number of simulations act to its disadvantage. Hence, a new methodology was proposed in this study by which the statistical properties of the dependent variables of the considered hydrologic problem may be obtained in only one single simulation.

The proposed FPE methodology derived in this study involved upscaling the governing stochastic differential equations by developing their corresponding Lagrangian–Eulerian Fokker–Planck Equation (LEFPE), thus transforming the original stochastic equations into the framework of a deterministic differential equation. The deterministic LEFPE that describes the time–space evolution of the probability density function of the unsteady open-channel flow state variables, was developed following the method in Kavvas (2003) after the governing Saint Venant equations were transformed into their characteristic form by using the method of characteristics. Through simplifications, this LEFPE was reduced to a classical FPE that could be solved deterministically for the evolution of the probability density of the state variables of the system. The obtained linear FPE was, then, discretized in an implicit manner following Chang and Cooper (1970). This provided the equations that may be solved numerically, through only one simulation, in order to determine the ensemble behavior and variability of a system described by the stochastic open-channel flow equations. The application and validation of this methodology, applied to an open-channel flow problem with an uncertain roughness coefficient, is provided in a companion paper by Dib and Kavvas (2018), in which the statistical results of the proposed FPE methodology are compared against the results obtained by the MC approach.

The open-channel flow problem considered in this study was for a
rectangular, prismatic channel under an uncertain roughness coefficient.
However, the proposed methodology can be expanded to problems which assume
uncertainties that arise from other flow or channel parameters. For
instance, when parameters such as the channel bed slope or the channel width
are assumed to be uncertain, their corresponding representations in the
equations will simply have to be included in the expectation and variance
expressions of the advection and diffusion coefficients for the

This study involved the theoretical development and derivation of a new proposed methodology for the stochastic solution of unsteady open-channel flow. No data were used in the derivation process.

The authors declare that they have no conflict of interest.