Ensemble modeling of stochastic unsteady open-channel flow in terms of its time-space evolutionary probability distribution : numerical application

The characteristic form of the Saint-Venant equations is solved in a stochastic setting by using a newly proposed Fokker–Planck Equation (FPE) methodology. This methodology computes the ensemble behavior and variability of the unsteady flow in open channels by directly solving for the flow variables’ time-space evolutionary probability distribution. 10 The new methodology is tested on a stochastic unsteady open-channel flow problem, with an uncertainty arising from the channel’s roughness coefficient. The computed statistical descriptions of the flow variables are compared to the results obtained through Monte Carlo (MC) simulations in order to evaluate the performance of the FPE methodology. The comparisons show that the proposed methodology can adequately predict the results of the considered stochastic flow problem, including the ensemble averages, variances, and probability density functions in time and space. Unlike the large 15 number of simulations performed by the MC approach, only one simulation is required by the FPE methodology. Moreover, the total computational time of the FPE methodology is smaller than that of the MC approach, which could prove to be a particularly crucial advantage in systems with a large number of uncertain parameters. As such, the results obtained in this study indicate that the proposed FPE methodology is a powerful and time-efficient approach for predicting the ensemble average and variance behavior, in both space and time, for an open-channel flow process under an uncertain roughness 20 coefficient.


Introduction
One of the most important types of unsteady open-channel flow problems is that of flood routing (Scharffenberg and Kavvas, 2011).This problem considers an initially uniform flow rate through an open channel, after which a flood wave enters upstream of the channel reach, and is translated downstream.The routing process, which predicts the spatial shape and temporal development of this flood wave as it traverses downstream (Viessman et al., 1977), involves solving for two dependent flow variables as a function of time and space through the river reach: velocity and depth, or discharge and depth.Solving for these two unknowns by means of the hydraulic routing technique involves using two governing equations (Chow, 1959), the continuity equation and the momentum equation, which are jointly known as the Saint-Venant equations (Sturm, 2001).
While it may be possible to assume that the flow and channel parameters are deterministic so as to obtain a deterministic solution to the Saint-Venant equations, such parameters usually exhibit high uncertainty in the real world (Gates and AlZahrani, 1996;Ercan and Kavvas, 2012a).In fact, uncertainties in these parameters may originate from several factors, including a channel's physical conditions and geometric parameters, its upstream boundary and initial conditions, and any lateral inflows (Chow, 1959;Sturm, 2001;Liang and Kavvas, 2008;Ercan and Kavvas, 2012a).As such, the parameters may be formulated as random functions, thus becoming spatially and/or temporally random, in which case the system can no longer be assumed deterministic.This necessitates the solution of the governing equations within a stochastic framework, from which a quantitative description Published by Copernicus Publications on behalf of the European Geosciences Union. of the ensemble behavior and variability of the process is obtained.In this manner, the two dependent flow variables are solved for their statistical properties, not for their deterministic values, at designated time-space positions.
Among the available approaches that can be used to solve the Saint-Venant equations within a stochastic framework, the Monte Carlo (MC) approach is one of the most wellknown due to its abundant use in simulating differential equations with stochastic parameters (Freeze, 1975;Smith and Freeze, 1979;Bellin et al., 1992).It is also generally accepted as the most robust approach for uncertainty evaluation, as well as the benchmark for comparing other new methods (Scharffenberg and Kavvas, 2011).However, one of the main disadvantages of the MC approach is its computational expense, which results from the large number of simulations that it usually involves.
In order to circumvent having to solve the Saint-Venant equations repeatedly for a large number of times, this study uses a new methodology that solves for the time-space evolutionary probability distribution of the unsteady open-channel flow process in only one simulation.From this probabilistic solution, one can then obtain the ensemble mean and ensemble variance of the process as they evolve in time and space.This new methodology, which is proposed, explained, and derived in the companion paper by Dib and Kavvas (2018), makes use of the ensemble averaging technique developed in Kavvas (2003) to obtain a Fokker-Planck Equation (FPE) that specifically describes an unsteady openchannel flow process.Some other hydrologic processes have been successfully simulated by following a similar procedure, which involved applying the ensemble averaging technique of Kavvas (2003) to their corresponding governing equations and obtaining the FPEs specific to their case.These include unsaturated water flow (Kim et al., 2005a), rootwater uptake (Kim et al., 2005b), 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).
Note that in addition to producing the statistical properties in a computationally efficient manner through one simulation, the FPE methodology developed for the unsteady open-channel flow process directly solves for, and is linear in, the probability density of the dependent variables.Moreover, while this methodology assumes a finite correlation time for the considered process, it does not make any linearization assumptions and it does not have limitations on the working range of the parameter space.
Therefore, in the wake of the preceding discussions, the first objective of this study is to use the FPE methodology derived in the companion paper for the unsteady open-channel flow process (Dib and Kavvas, 2018) and to apply it to a representative stochastic unsteady open-channel flow problem in order to solve for the probability density of the state variables of the flow process and to provide a quantitative de-scription of the expected behavior and variability of the system in one single simulation.The second objective is to evaluate the performance of the proposed methodology and to validate its results by comparing the statistical properties of the flow variables computed by the FPE methodology against those calculated by the MC approach.
2 Saint-Venant equations: characteristic form and ensemble-averaged form The Saint-Venant equations solved in this study are written for the unsteady open-channel flow of an incompressible fluid in a rectangular, prismatic channel with no lateral inflow.The method of characteristics is used to transform these equations from two partial differential equations into a system of four ordinary differential equations (ODEs).These four ODEs include two characteristic equations describing the two characteristic paths, and two compatibility equations which must be satisfied along their corresponding characteristic path.The characteristic form of the Saint-Venant equations is shown below (Sturm, 2001): Flow process condition to be satisfied along Flow process condition to be satisfied along C where V is the average flow velocity, c is the wave celerity which is equal to √ gy for a rectangular channel where y is the flow depth, x is the position, t is the time, S 0 is the slope of the channel bottom, S f is the friction slope, and g is the acceleration of gravity.Note that S 0,i denotes S 0 (x i , t); the same applies to S f .The positive and negative characteristic curves are defined as C 1 and C 2 , respectively, and the variables or derivatives corresponding to C 1 and C 2 are denoted by the subscripts 1 and 2, respectively.
The Lagrangian-Eulerian form of the Fokker-Planck equation (LEFPE), corresponding to the Saint-Venant equations, that can solve for the multivariate probability density function (PDF) of the hydrologic state variables of an unsteady open-channel flow problem was obtained through the following steps: (i) using the two substitutions V + 2c = α and V − 2c = β on Eqs.(1) to ( 4), (ii) applying the technique in Kavvas (2003) while assuming the uncertainty arising from the Manning's roughness coefficient, and (iii) per-forming several simplifying assumptions.The result is shown in Eq. ( 5), while its detailed derivation can be found in the companion paper by Dib and Kavvas (2018). (5) In Eq. ( 5), b denotes the width of the channel, n denotes Manning's roughness coefficient, and k denotes the conversion factor between SI and US units for Manning's formula.The remaining variables are as defined previously.Note that in Eq. ( 5), and other later equations, P (x 1 , x 2 , α, β, t) is sometimes substituted by P for simplicity.Moreover, note that Eq. ( 5) resembles an advection-diffusion equation, in which the first four bracketed terms multiplied by P on their left-hand sides resemble the advection coefficients, and the last four bracketed terms resemble the diffusion coefficients.Hence, after denoting the advection coefficients by F and the diffusion coefficients by D, Eq. ( 5) was simplified into Eq.( 6) below, which is the final analytical form of the proposed FPE methodology.
An appropriate numerical scheme to solve Eq. ( 6) was derived following Chang and Cooper (1970), as discussed in detail in the companion paper.This scheme involves implicitly discretizing the equation, while noting that the dependent variables are to be solved at the intersection of the characteristic curves C 1 and C 2 , which renders As such, the discretized version of Eq. ( 6) was derived in the companion paper to be as follows: where h: 0, 1, 2, . . ., N H denotes the domain of x; k: 0, 1, 2, . . ., N K denotes the domain of α; l: 0, 1, 2, . . ., N L denotes the domain of β; and n: 0, 1, 2, . . .denotes the domain of time t.The expression for the λ α parameter is given in Eq. ( 8) below, while the expressions for the other λ parameters can be written in a similar manner.
A. Dib and M. L. Kavvas: Ensemble modeling of stochastic unsteady open-channel flow -Part 2 where i: 0, 1, 2 . . ., N I denotes the domain of x 1 in the direction of the C 1 curve; j : 0, 1, 2 . . ., N J denotes the domain of x 2 in the direction of the C 2 curve.Equation ( 7) can be solved implicitly to compute the joint PDF of the state variables within the x-α-β domain.From this, a quantitative probabilistic description of the ensemble behavior and variability can be determined for the hydrologic system represented by the stochastic Saint-Venant equations with uncertain parameters.
3 Application of the Monte Carlo approach and the new Fokker-Planck equation methodology to a hydraulic routing problem A hypothetical hydraulic routing problem will be the illustrative example which the new proposed methodology will be tested on.The hypothetical problem involves a river reach that is 2.70 km long, sloping at 0.0015 throughout the whole reach, having a rectangular cross section with a constant width of 6.1 m, and having no lateral inflow.Initially (at t = 0), the river is assumed to have a steady uniform flow of 15.5 m 3 s −1 throughout the reach.At t > 0, the upstream flow increases linearly to reach 56 m 3 s −1 at t = 20 min, then decreases linearly back to reach the initial flow value at t = 60 min, after which it remains at the same constant flow value for t > 60 min.For the purpose of this study, the problem is made stochastic through the uncertainty in the Manning's roughness coefficient, whose statistical characteristics must be estimated.All simulations for the MC approach and the FPE methodology were run on a computer having 16 GB of RAM and an Intel i7 processor with four cores, each core having a base frequency of 2.40 GHz and a maximum frequency of 3.40 GHz.

Estimating the statistical characteristics of Manning's roughness coefficient
Manning's roughness coefficient (n) represents the resistance of the bed of a channel to the flow of water in it (Chow, 1959).Generally, a higher n value defines a greater resistance to the flow.The value of this coefficient depends on several factors, including surface roughness, vegetation, obstructions, and channel irregularity and alignment (Chow, 1959), all of which may exhibit considerable variability along the length of an open channel.It is, therefore, crucial to account for such variability in order to better represent the behavior of the unsteady open-channel flow system being solved.
The published and technical studies with sizeable datasets to address the variability of n in such open channels are few.However, some studies have provided median and range values, while others have attempted to fit different probability distributions to the data (Gates and AlZahrani, 1996).Using such information, Manning's n was assumed for this study to be a random variable with a normal distribution having a mean of 0.035 and a standard deviation of 0.005.The mean and standard deviation were chosen in a way that provides realistic values and distributions of n that fall within the ranges and statistics provided by Gates and AlZahrani (1996) and that conform to the table of typical values in Chow (1959).Moreover, the selected mean and standard deviation allowed the generation of n values which never fell below 0.01, thus complying with the fact that the roughness coefficients for flows in natural streams and excavated channels are always greater than 0.01 (Chow, 1959).As such, no truncation or discarding of any generated n values was required.The chosen PDF for the roughness coefficient was used by the MC approach as well as by the new FPE methodology when solving for the ensemble behavior and variability of the hypothetical flow problem.

Application of the Monte Carlo approach
For the hypothetical routing problem of this study, the MC approach requires repeatedly solving the Saint-Venant equations in a deterministic manner for a large number of different roughness coefficient (n) realizations.To deterministically solve the Saint-Venant equations in their full form, several numerical techniques have been developed because their analytical solution has not been possible due to the presence of nonlinear terms (Sturm, 2001;Chaudhry, 2008).For this study, the characteristic form of the Saint-Venant equations was discretized in an explicit manner by substituting the time derivatives with their first-order finite-difference forms, as detailed in several references, e.g., Viessman et al. (1977) and Sturm (2001).The values of the dependent variables at the new time steps were computed at the points of intersection of the positive and negative characteristic curves, which rendered the final solution on an irregular x-t grid.The solution was then interpolated onto a rectangular grid, with a x of 75 m and a t of 3 min.The simulations were performed using a parallelized process which optimized the computational time by running the simulations over the total number of available cores (with no hyper-threading).
The stability of the numerical method being used was ensured by checking the Courant condition (Sturm, 2001).Furthermore, two boundary conditions, one at each end of the reach, and two initial conditions were defined in this study, since the problem deals with the subcritical unsteady nonuniform flow case (Sturm, 2001).As initial conditions, the discharge at every location along the river was provided (taken as 15.5 m 3 s −1 , as explained in the problem description), and the flow was assumed to be initially uniform and steady.As for the boundary conditions, the flow hydrograph at the channel entrance was given, while at the downstream end, the channel was assumed to be hydraulically long so that the flow can be taken as normal flow, thus satisfying Manning's equation.As such, the downstream boundary condition was chosen as the depth-discharge relationship represented by Manning's equation.This equation, along with the equations corresponding to the positive characteristic (i.e., Eqs. 1 and 2), was used to compute the flow variables at the downstream boundary.
Following the preceding discussion, the Saint-Venant equations were deterministically solved for a total of 1000 times, each time using a different realization of n that was generated based on the PDF chosen in Sect.3.1.While a lower number of realizations may have been sufficient for accurate computations of the first moment of the flow variables, it would not have been sufficient for the accurate computations of the second moment.In fact, the standard deviation of the flow discharge computed using 50, 100, 200, and 500 realizations showed absolute relative differences that reached 65, 35, 20, and 15 %, respectively, when compared to the results of the 1000-realization run.Therefore, the number of realizations in this study was selected to be large enough to numerically approximate, with sufficient accuracy, both the first and the second moments of the stochastic quantities of the problem at hand.

Application of the proposed Fokker-Planck equation methodology
Unlike the MC approach, the FPE methodology aims at solving for the ensemble behavior and variability of the stochastic unsteady open-channel flow system in one shot by computing the PDF of its dependent variables over time and space.
Solving the hypothetical routing problem following the FPE methodology involved solving Eq. ( 7).Since this equation is the result of implicit discretization, it is unconditionally stable and requires no constraint on the size of the time step for its stability.However, to obtain sufficiently accurate solutions, the time step must still be limited with the Courant condition (Zheng and Bennett, 2002).Therefore, at every time position, the multidimensional Courant condition was checked in order to determine the appropriate size of the next time step.Furthermore, to correctly solve the hydraulic routing problem, both the initial and boundary conditions of the problem in the physical space must be accurately represented in the probability space for use in the FPE methodology.Note that when the initial condition in the physical space is taken as deterministic, the Dirac delta function (δ(s)) can be used to represent it in the probability space.Assuming that H is a vector of m state variables (m dimensional), and H 0 is the corresponding vector of initial conditions, the initial condition PDF of H can be written as For the purpose of this study, H corresponds to a point in the x-α-β domain (3 dimensional), and H 0 represents the initial condition that defines the deterministic flow at time t = 0 in the x-α-β domain.However, since it is not possible to numerically represent the Dirac delta function as a function with infinite value at only one position, it was estimated in the probability domain as a function with a very high value over a small bottom width, while preserving an area of unity.
Concerning the boundary conditions, they are usually divided for FPEs into two categories: accessible and inaccessible.Inaccessible boundaries are defined as those boundaries that could never be reached if the process starts from any interior point of the domain.On the other hand, for accessible boundaries, there is a positive probability that these boundaries will be reached from the interior of the domain within a finite amount of time (Feller, 1954).Accessible boundaries can be further subdivided into absorbing and reflecting, or no-flux, boundaries.
Note that since the FPE can be considered as the conservation equation for probability mass, a probability mass of unity needs to be conserved in the probability domain of the system.Therefore, using a reflecting (no-flux) boundary condition would be the most suitable choice for this study in order to ensure the completeness of the probability domain and to prevent any probability mass from leaving the domain (Gardiner, 1985).Such a condition was used to describe the boundaries of both the α and β dimensions, noting that the minimum and maximum boundaries of both α and β were chosen to be far enough so that they would encompass all possibilities that could occur for the considered routing problem.As for the upstream boundary, recall that the discharge hydrograph was assumed to be known upstream, in which case the probability densities at the upstream boundary of the x dimension were known for all t > 0. Finally, the downstream boundary in the x dimension was formulated to replicate that of the MC model, and it was extended downstream much further than the required length of the reach so as to eliminate any of its effect on the numerical solution.

Numerical results and discussion
Both the MC approach and the proposed FPE methodology were used to solve the stochastic Saint-Venant equations of the hypothetical hydraulic routing problem presented in Sect.3, in order to determine the ensemble behavior of the system and the statistical distributions of the flow variables.While the state variables directly solved for were the velocity and depth (or celerity), the discharge was easily computed from these two variables.
A plot of the ensemble average discharge over time and space computed by the FPE methodology can be seen in  Fig. 1 alongside a plot of the same results that were obtained by the MC simulations.From this figure, it is clear that the ensemble average discharge computed by the FPE methodology resembles the one obtained from the MC simulations quite well, while showing the same behavior and evolution of the mean discharge in both time and space as a result of the applied upstream wave.From the average behavior of the system, both plots show that the wave that was initiated upstream is observed to be transmitted downstream through the reach.In this case, the discharge at every location is seen to increase from the initial value of 15.5 m 3 s −1 to some peak value, after which it decreases back again to the initial flow value.However, it is noticeable in both plots that the aver-age peak discharge becomes lower at locations further downstream.This decrease comes as a result of the dissipation of energy through viscous effects as the water flows along the reach (Jeppson, 2011).In fact, the shearing stresses due to the vertical velocity gradients within the water lead to a loss of the fluid energy (kinetic and potential) as non-recoverable energy (e.g., increase in temperature), causing the peak velocity and discharge to decrease along the reach.
For a clearer comparison between the FPE and MC mean discharge results, cross sections of both plots from Fig. 1 at specific times (Fig. 2) and at specific channel locations (Fig. 3) were compared individually.From Fig. 2, it is clear that the change in the ensemble average discharge as a func-Hydrol.Earth Syst.Sci., 22,[2007][2008][2009][2010][2011][2012][2013][2014][2015][2016][2017][2018][2019][2020][2021]2018 www.hydrol-earth-syst-sci.net/22/2007/2018/ tion of position (location within the channel reach along the flow direction) as computed by the FPE methodology almost perfectly replicates the corresponding results obtained by the MC simulations.In both methods, the effect of the wave is clear as it causes an increase in the discharge first upstream, and then throughout the reach over time.Comparing the mean discharge at specific positions over time (Fig. 3) reveals that the FPE methodology also predicts the temporal change in the mean discharge very well at different channel locations.While a slight decrease in the discharge can be seen at early times of Fig. 3c and d, which may be attributed to the numerical discretization, the timing and positions of the peak discharges predicted by the FPE methodology are very similar to those of the MC approach in all plots of Fig. 3, with a maximum relative difference of only around 6 %.Similar results were also plotted for the ensemble average flow depth.Figure 4 shows the mean river channel flow profile at different times as computed by the FPE and the MC methods.This figure shows how the applied upstream wave affects the depth profile of the river at different times, which is similarly predicted in both the FPE and MC simulations.Compared to the MC results, the ensemble average depth is predicted quite well by the FPE methodology, which only slightly overestimates the mean depth values at times between t = 30 min and t = 45 min.A comparison of the changes in depth as a function of time (Fig. 5) reveals that the FPE methodology provides a good match to the MC results while showing the same peaking pattern as the wave reaches a specific channel location.The timings of the peaks of the FPE methodology greatly match those of the MC simulations.Similarly to Fig. 4, a slight overestimation can be noticed from the FPE methodology, especially around the peak depths, but with a maximum relative difference of only around 7.5 %.
As for the velocity, Figs. 6 and 7 show the same two kinds of plots as before by comparing the mean ensemble average velocity results of the FPE methodology to those of the MC approach.Figure 6 shows that the mean velocity computed by the FPE methodology provides a good match to the MC results, with a slight underestimation.The pattern of change in the mean velocity as a function of position that is presented by the FPE methodology is very similar to that shown by the MC simulations, which demonstrates the effects of the applied upstream wave as it travels downstream with time.Figure 7 also reveals a good match between the two modeling approaches, with the FPE methodology providing the same pattern of velocity change over time as computed by the MC simulations.Similarly, a slight underestimation of the FPE mean velocity is also seen in Fig. 7, but the maximum absolute relative difference between both methods was only around 11 %.Therefore, Figs. 1 to 7 show the capability of the FPE methodology in predicting the ensemble averages of the flow variables quite well over space and time, which in turn provides support and validation to the simplifications and assumptions applied to the FPE methodology.In a similar manner to the ensemble averages, the relative performance of the FPE methodology in predicting the system's variability was examined, this time by checking the standard deviations.Figure 8 shows a comparison of the standard deviation of the flow discharge over space and time as computed by the FPE methodology and by the MC simulations.Both plots of this figure reveal that the standard deviation experiences two triangular areas of high values, the earlier in time being generally higher than the later, and both areas showing a stronger intensity further downstream.While the general resemblance of the FPE plot to the MC plot is good, the second area of elevated standard deviations is more compact in the FPE results and shows slightly lower values than those of the MC results.In an attempt to study such differences more closely, cross sections of both plots from Fig. 8 at specific times and at specific channel locations were compared individually (Figs. 9 and 10), in a manner similar to the ensemble average plots discussed previously.
Looking at the change in the standard deviation of the discharge over position, Fig. 9 reveals how the movement of the wave in the downstream direction causes an increase in the standard deviation of the discharge, thus increasing the variability in the discharge.After the whole channel is affected by the moving wave, the variability is seen to be highest at the downstream end of the channel (e.g., Fig. 9b-d).As the wave subsides, the standard deviation starts to decrease (Fig. 9c and d).Comparing the results of the FPE methodology to those of the MC approach shows that the FPE results correctly predict the patterns of change in the standard deviation as a function of position at each of the different times (Fig. 9).While the values of the standard deviations are slightly overestimated in some of the FPE results, especially at t = 30 min, the FPE results are still quite close to the MC results.
Concerning the changes in the standard deviation of the discharge over time, Fig. 10 shows that the standard deviation forms an M-shaped pattern at all the plotted locations, which can also be deduced from Fig. 8.For the MC results, the magnitudes of the standard deviation forming this M-shaped pattern are observed to be generally increasing as a function of the distance away from the upstream boundary.It is noticeable that the results of the FPE methodology also predict an M-shaped pattern and an increase in the magnitude of the standard deviation of the discharge as a function of location.However, the standard deviation of the FPE methodology shows an offset at time t = 0 when compared to the MC simulations, which can also be noticed at the upstream positions of Fig. 9. Recall that the initial and upstream flow discharge is assumed to be known.Nonetheless, a single known value of the flow discharge, when joined with a spread of roughness coefficients due to the uncertainty involved, leads to an unavoidable spread in the velocity and depth values.Since α and β are functions of the velocity and celerity (and in turn, depth), this spread is translated onto the α-β plane of the FPE methodology.As a result, with an uncertain rough-ness coefficient, the only way to numerically represent a deterministic discharge in the α-β plane is to have a spread of probability mass over the values involved.The existence of this spread on a non-continuous, discretized α-β plane may have had the most contribution to the offset of the standard deviation at the initial times and positions of Figs. 9 and 10.
Nonetheless, when looking at the standard deviation after the initial times, one can see that the timings, shapes, and magnitudes of the first peak for the MC results are generally matched well by the FPE results as shown by the plots of Fig. 10.Moreover, the locations and magnitudes of the second peak are well represented by the FPE results in Fig. 10a  and b.However, differences arise in the position and magnitude of the second peak for locations further downstream (Fig. 10c and d).In fact, while moving further downstream, the second peak is observed to be underestimated by the FPE results and slightly shifted back in time when compared to the MC simulations.Moreover, the standard deviation following the second peak shows a faster decline in time for the FPE methodology than for the MC approach, while the local minimum between the two peaks is seen to be overestimated by the FPE results.
Several factors may have caused the discrepancies in the standard deviation of the discharge as computed by the FPE methodology.Among those factors may be the approximations and assumptions that were applied in order to simplify the FPE methodology to an easily solvable degree.These approximations have simplified the computations of the drift and diffusion coefficients of the equation, thus causing some discrepancies in the movement and the spread of the probability mass within the domain.Other factors may include the numerical method selected for computing the FPE methodology, as well as the associated spatial and temporal discretizations that were used.All of these may have had an effect on the representation of the probability flow within the domain, thus leading to some discrepancies in the variability results.Such discrepancies may be faced when using the FPE methodology in engineering flow applications, such as flood forecasting and flood control.The variability of the flow in flood forecasting applications, for example, may be underestimated at the downstream end of the reach, specifically during the later time periods.This would impact the range of flows that are forecasted to occur at the downstream end.
Nonetheless, it may still be inferred that the FPE methodology performs satisfactorily in predicting the variability of the discharge in both space and time.In fact, the problem being solved is highly nonlinear, and involves a large uncertainty in one of its parameters.As such, estimations in the mean and variance are even more difficult to accurately quantify.The FPE methodology was capable of not only correctly matching the general patterns of the spatial and temporal changes in the discharge standard deviation but also correctly providing the standard deviation values within a range that is very similar to the range computed by the MC simulations.
Concerning the standard deviations of the velocity and depth, it is important to note that their behavior over position is somewhat different from that of the flow discharge.In fact, the standard deviations of the velocity at the same four time positions of Fig. 9 seem to be relatively constant at each time position, having a value between 0.015 and 0.02 m s −1 .On the other hand, the standard deviation of the depth showed a greater range of values at each time position, as a function of location, with values ranging between 0.15 and 0.5 m.When looking at the standard deviations as a function of time at the same four locations of Fig. 10, both standard deviations seem to show that their values increase to reach a maximum and then decrease to levels similar to original levels, not unlike their corresponding ensemble average plots over time.Again, the range of change in the standard deviation of the velocity is much smaller (0.015 to 0.02 m s −1 ) than that of the depth (0.15 to 0.5 m).Note that the relative differences in the FPE results when compared to the MC results reach up to 23 and 29 % for velocity and depth, respectively.
A final comparison between the FPE methodology and the MC approach is to compare the probability density functions (PDFs) of the flow discharge at different times and locations.The large number of simulations (1000) done using the MC approach provided an equal number of flow discharge results at each specific x-t position in the space-time plane.These values were then used in order to estimate the PDF of the discharge at that specific x-t position.The FPE methodology, on the other hand, directly solved for the evolution of the PDFs of the state variables through time and space.
However, recall that the FPE was solved in the x-α-β domain.With the discharge being a function of both α and β, the PDF of the discharge at each specific x-t position was deduced from the x-α-β PDF provided by the FPE methodology.
For comparison, the PDFs of the flow discharge obtained from the FPE and MC methods were plotted at three different channel locations (x = 900, 2250, and 2700 m) and at four different times for each location (t = 15, 30, 45, and 60 min), as shown in Fig. 11.For most of the plots, the peak values of the PDFs predicted by the FPE methodology are similar to those computed by the MC approach.However, taking Fig. 11g as an example, while the peak of the MC PDF has a value of around 0.35, that of the FPE methodology has a value of around 0.065.This is because the PDF of the FPE methodology has a larger spread, which causes the reduction in its PDF values to preserve an area of unity.It may be noted that this plot corresponds to a point close to the local minimum of the plot in Fig. 10b, in which the FPE methodology is seen to predict a larger standard deviation than the MC approach.Therefore, while both of the PDFs, corresponding to the FPE and MC approaches, seem to be located within the same range and thus providing similar expected values of the discharge, the greater variability of the FPE PDF causes its lower peak values.
Many different factors may have affected the evolution and calculation of the PDFs for the flow discharge in space and time.Such factors may even include the step sizes in the α and β directions whose values directly affect the computation of the discharge PDFs.However, the PDFs predicted by the FPE methodology are generally seen to be following similar trends to those computed by the MC approach (Fig. 11) while also satisfactorily predicting the ranges and locations of the PDFs (along the x-axis).Thus, with all the variability encompassing the routing problem considered, and with all the assumptions and approximations used during the application of the FPE methodology, the probability densities predicted by the proposed FPE methodology are considered to be rather encouraging.
It should be noted that these FPE results required significantly less time for computation as opposed to the MC results.Recall that the 1000 MC simulations were parallelized and run over all four cores (with no hyper-threading), thus noticeably reducing the computational time as compared to an un-parallelized run.With such parallelization, the MC simulations ran for over 2 days.On the other hand, the results of the FPE methodology, which was not parallelized, were obtained in about 7 h.
If we observe the computational times of the implicit numerical solution of the FPE methodology, the portion of the simulation requiring the greatest time is filling out the coefficient matrix, especially for small α and β discretizations.Parallelizing this portion over the four cores would allow one to considerably reduce the time to fill out the coefficient matrix, thus reducing the total computational time of this method.Without the parallelization of the FPE methodology, its one simulation may still not seem to provide an immense advantage when only one uncertain parameter is involved, especially with the possibility of parallelizing the MC simulations among a much larger number of cores.Nonetheless, when the problem being solved involves a greater number of uncertain parameters and boundary conditions, or even a larger system, such an advantage may prove to be crucial.In fact, the computational expense of the MC simulations for such a case would be expected to increase exponentially due to the higher number of simulations needed to maintain the desired accuracy in the results, thus significantly increasing the computational time regardless of parallelization.On the other hand, such additional uncertainties can be easily imple-  mented into the FPE methodology by making simple changes and additions that will be reflected in Eq. ( 5), after which the FPE would be solved following similar steps as discussed for this study, with minimal implications on the computational expense.

Summary and conclusions
This study applied the proposed FPE methodology derived in the companion paper by Dib and Kavvas (2018) to a stochastic unsteady open-channel flow problem, with an uncertain roughness coefficient.The equations used to describe the open-channel flow problem were the Saint-Venant equations, transformed into their characteristic form by the method of characteristics.The proposed FPE methodology was applied in order to solve for the probability density of the flow state variables (velocity and depth/celerity, as well as discharge) and to provide a quantitative description of the expected behavior and variability of the stochastic system in one single simulation, as opposed to the large number of simulations usually performed by the MC approach.The unsteady open-channel flow problem was also solved using the MC approach in order to use its results for evaluating the performance of the proposed FPE methodology.
Comparisons of the FPE results to those of the MC simulations revealed the effectiveness of the proposed FPE methodology in describing the ensemble behavior and variability of the stochastic flow problem.In fact, the FPE methodology was found to replicate the ensemble average discharge of the MC simulations quite well in both space and time.In addition, it was also capable of effectively representing well the temporal and spatial change in the ensemble average depth and the ensemble average velocity.Furthermore, this method provided a good representation of the patterns and ranges of the standard deviation of the discharge over time and space, and showed a satisfactory prediction of the spatial and temporal trends and ranges of the flow discharge PDFs.These encouraging results were obtained despite simplifications applied to the FPE methodology.Furthermore, with the FPE approach, the computational time was significantly less than the time taken by the MC approach, and the FPE methodology results were obtained by running only one simulation, as opposed to the large number of simulations performed by the MC approach.Such an advantage becomes even more prominent with a greater number of uncertain parameters and boundary conditions, in which case the computational expense of the MC simulations that is needed to preserve the desired accuracy would exponentially increase.On the other hand, only simple adjustments would be required for the FPE, which could then be solved as was done in this study, with minor implications on its computational expense.Therefore, the results obtained in this study indicate that the proposed FPE methodology may be a powerful and timeefficient approach for predicting the ensemble average and variance behavior, in both space and time, for the unsteady open-channel flow process under an uncertain roughness coefficient, hence being an approach that would be essential for engineering flow problems.
While the FPE methodology satisfactorily described the ensemble average and variability of the open-channel flow system in this study, this methodology is open to improvements especially with regard to reducing any discrepancies in its numerical results.Running a more comprehensive version of the FPE methodology, by including only some of the simplifying assumptions used in this study, may be one option.Another option may involve using a higher-order and more accurate numerical scheme for the discretization of the multidimensional FPE.As such, numerous opportunities present themselves for future research within this topic, all of which would be of great benefit in the further improvement of the proposed methodology.

Figure 1 .
Figure 1.Ensemble average discharge over channel position and time obtained by (a) the FPE methodology and (b) the MC approach.

Figure 2 .
Figure 2. Comparison of the ensemble average discharge obtained by the FPE methodology and the MC approach as a function of channel location (Position), at different times.

Figure 3 .
Figure 3.Comparison of the ensemble average discharge obtained by the FPE methodology and the MC approach as a function of time, at different channel locations.

Figure 4 .
Figure 4. Comparison of the ensemble average flow depth obtained by the FPE methodology and the MC approach as a function of channel location, at different times.

Figure 5 .Figure 6 .
Figure 5.Comparison of the ensemble average flow depth obtained by the FPE methodology and the MC approach as a function of time, at different channel locations.

Figure 7 .Figure 8 .
Figure 7.Comparison of the ensemble average velocity obtained by the FPE methodology and the MC approach as a function of time, at different channel locations.

Figure 9 .
Figure 9.Comparison of the standard deviation of the flow discharge obtained by the FPE methodology and the MC approach as a function of channel location, at different times.

Figure 10 .
Figure 10.Comparison of the standard deviation of the flow discharge obtained by the FPE methodology and the MC approach as a function of time, at different channel locations.

Figure 11 .
Figure 11.Comparison of the probability density functions of the flow discharge obtained by the FPE methodology and the MC approach, plotted at different times and channel locations.