An Efficient Solver for Cumulative Density Function-based Solutions of Uncertain Kinematic Wave Models

TL;DR

This paper introduces a computationally efficient numerical framework using the CDF method and Monte Carlo simulations to determine the probability distribution of uncertain kinematic wave models, validated against the Saint-Venant equation.

Contribution

The paper presents a novel, efficient solver for the CDF-based approach to uncertain kinematic wave models using the method of characteristics and Monte Carlo simulations.

Findings

The proposed method is significantly faster than direct Monte Carlo simulations.

The solver accurately reproduces the probability distribution of the system state.

Validation against the Saint-Venant equation confirms robustness and accuracy.

Abstract

We develop a numerical framework to implement the cumulative density function (CDF) method for obtaining the probability distribution of the system state described by a kinematic wave model. The approach relies on Monte Carlo Simulations (MCS) of the fine-grained CDF equation of system state, as derived by the CDF method. This fine-grained CDF equation is solved via the method of characteristics. Each method of characteristics solution is far more computationally efficient than the direct solution of the kinematic wave model, and the MCS estimator of the CDF converges relatively quickly. We verify the accuracy and robustness of our procedure via comparison with direct MCS of a particular kinematic wave system, the Saint-Venant equation.

Figures

Tables

Table 1. Table 1: l 2 subscript 𝑙 2 l_{2} error between the numerical result and the exact solution ( 13 ) at time t = 0.1 𝑡 0.1 t=0.1 .

$\Delta t$	$ϵ$	Convergence Rate
$0.1$	$8.95\times {10}^{-3}$	—
$0.05$	$1.11\times {10}^{-3}$	$3.02$
$0.025$	$1.37\times {10}^{-4}$	$3.01$
$0.0125$	$1.71\times {10}^{-5}$	$3.00$

Table 2. Table 2: The CPU time of the direct simulation (MCS) and the fine-grained CDF approach (CDF) for a single realization, using a workstation with an Intel (R) Core (TM) i5-4210H 2.90 GHz with 8.00 GB RAM.

	$S=0$	$S=1$	$S=x$
MCS (seconds)	46.22	45.05	45.99
CDF (seconds)	2.00	1.98	2.00

Table 3. Table 3: Comparison between the WENO-Roe numerical result and the exact solution ( 13 ) at time t = 0.1 𝑡 0.1 t=0.1 .

$\Delta x$	$ϵ$	Convergence Rate
0.05	$2.22\times {10}^{-5}$	None
0.025	$1.02\times {10}^{-6}$	4.44
0.0125	$4.59\times {10}^{-8}$	4.47
0.00625	$1.97\times {10}^{-9}$	4.54

Equations

$$\mathbf{q}=\mathbf{q}(k,\mathbf{z}),$$

$$\frac{\partial k}{\partial t}+\nabla \cdot \mathbf{q}=S,\mathbf{x}\equiv (x_1,x_2,x_3)\in \mathcal{D},t>0,$$

$$k_{\mathrm{in}}(\mathbf{x})=k(\mathbf{x},t=0),$$

$$k_{\mathrm{bx}}(t)=k(\mathbf{x}=\partial \mathcal{D},t),$$

$$\Pi (K;\mathbf{x},t)=\mathcal{H}[K-k(\mathbf{x},t)],K\in {\mathbb{R}}^{+},$$

$$\displaystyle\mathcal{H}(y)=\left\{\begin{array}[]{rl}0,&y<0\\ 1,&y\geq 0.\end{array}\right.$$

$$\mathbb{E}\Pi (K;\mathbf{x},t):={\int }_{-\infty }^{\infty }\mathcal{H}(K-k^{\prime })f_{\mathrm{k}}(k^{\prime };\mathbf{x},t)\mathrm{d}{k}^{\prime }=F_{\mathrm{k}}(K;\mathbf{x},t).$$

$$\frac{\partial \Pi }{\partial t}+\mathbf{v}(K;\mathbf{x},t)\cdot {\nabla }_{\mathbf{c}}\Pi =0,$$

$${\nabla }_{\mathbf{c}}={\left(\frac{\partial }{\partial x_1},\frac{\partial }{\partial x_2},\frac{\partial }{\partial x_3},\frac{\partial }{\partial K}\right)}^{\mathrm{T}},$$

$$\mathbf{v}(K;\mathbf{x},t)={\left(\frac{\partial q_1}{\partial K},\frac{\partial q_2}{\partial K},\frac{\partial q_3}{\partial K},\sum _{i,j=1}^{3,\mathrm{N}}\frac{\partial q_i}{\partial z_j}\frac{\partial z_j}{\partial x_i}+S\right)}^{\mathrm{T}},$$

$${\Pi }_{\mathrm{in}}=\Pi (K;\mathbf{x},0)=\mathcal{H}\left[K-k_{\mathrm{in}}(\mathbf{x})\right],$$

$${\Pi }_{\mathrm{bx}}=\Pi (K;\mathbf{x},t)=\mathcal{H}\left[K-k_{\mathrm{bx}}(t)\right],$$

$${\Pi }_{\mathrm{bk}}=\Pi (0;\mathbf{x},t)=0.$$

$$\frac{\mathrm{d}\mathbf{c}}{\mathrm{d}t}=\mathbf{v},\mathbf{c}=(x_1,x_2,x_3,K{)}^{\mathrm{T}},\mathbf{c}(t=0)={\mathbf{c}}_0.$$

$$\frac{\mathrm{d}\Pi }{\mathrm{d}t}=0.$$

$$F_{\mathrm{k}}(K;\mathbf{x},t)$$

$$\frac{\partial F_{\mathrm{k}}}{\partial t}+{\mathbf{v}}_{\mathrm{eff}}\cdot {\nabla }_{\mathbf{c}}{F}_{\mathrm{k}}={\nabla }_{\mathbf{c}}\cdot (\mathbf{D}{\nabla }_{\mathbf{c}}{F}_{\mathrm{k}}),$$

$$\frac{\partial k}{\partial t}+\frac{\partial q}{\partial x}=S,q=k^{\frac{1}{2}},x\in [0,2],t\in [0,+\infty )$$

$$\frac{\partial \Pi }{\partial t}+\frac{1}{2\sqrt{K}}\frac{\partial \Pi }{\partial x}+S\frac{\partial \Pi }{\partial K}=0.$$

$$\begin{array}{rl}& S(x,t)=2\pi \left[\sin (\pi (x+t))+1.1\right]\cos \left(\pi (x+t)\right)+\pi \left[\cos \left(\pi (x+t)\right)\right],\\ & k_{\mathrm{bx}}={\left[\sin \pi t+1.1\right]}^2,k_{\mathrm{in}}={\left[\sin \pi x+1.1\right]}^2,\end{array}$$

$$k_{\mathrm{exact}}(x,t)={\left[\sin \left(\pi (x+t)\right)+1.1\right]}^2.$$

$$ϵ:=\frac{1}{N}\sum _{i=1}^N{\left[k_{\mathrm{exact}}(x_i,t)-k(x_i,t)\right]}^2,$$

$$\begin{array}{rl}& S(x,t)=2z^2\pi \left[\sin \pi (x+t)+5\right]\cos \pi (x+t)+z\pi \cos \pi (x+t),\\ & k_{\mathrm{bx}}={\left[\sin \pi t+5\right]}^2,k_{\mathrm{in}}={\left[\sin \pi x+5\right]}^2,\end{array}$$

$$f_{\mathrm{z}}(Z;\mu ,{\sigma }^2)=\frac{1}{Z\sigma \sqrt{2\pi }}{\mathrm{e}}^{-{\left(\ln Z-\mu \right)}^2/2{\sigma }^2},$$

$$k(x,t,z)={\left[z\sin \pi (x+t)+5\right]}^2,$$

$$F_{\mathrm{k}}(K;x,t)$$

$$ϵ(x,t):=\frac{{\sum }_{i=1}^N{\left[F_{\mathrm{k},M}(K_i;x,t)-F_{\mathrm{k}}(K;x,t)\right]}^2}{{\sum }_{i=1}^N{F}_{\mathrm{k}}^2(K;x,t)},$$

$$\frac{\partial k}{\partial t}+\frac{\partial k}{\partial x_1}+\frac{\partial k}{\partial x_2}+\frac{\partial k}{\partial x_3}=\frac{\pi z}{3}[t\cos (\pi t{x}_1)+t\cos (\pi t{x}_2)$$

$$+t\cos (\pi t{x}_3)+x_1\cos (\pi t{x}_1)+x_2\cos (\pi t{x}_2)+x_3\cos (\pi t{x}_3)],$$

$$k=\frac{z}{3}\left[\sin (\pi t{x}_1)+\sin (\pi t{x}_2)+\sin (\pi t{x}_3)\right].$$

$$\frac{\partial k_1}{\partial t}+\frac{\partial k_2}{\partial x}$$

$$\frac{\partial k_2}{\partial t}+\frac{\partial k_1}{\partial x}$$

$$k_1(x,0)$$

$$k_2(x,0)$$

Full text

1 Introduction

Kinematic wave models (KW) [14, 15] provide an important mathematical tool to describe environmental flows, such as overland flow and erosion [24, 25]. Due to multiscale heterogeneity and insufficient site characterization, input parameters or fields into KW models often exhibit a high degree of uncertainty and are commonly modeled as random quantities. Such a probabilistic approach renders an otherwise deterministic kinematic wave equations as stochastic, and results in complete solutions that are in the form of probabilistic density functions (PDFs) or cumulative density functions (CDFs) of system states. This amounts to forward propagation of parametric uncertainty through the modeling process.

To quantify such predictive uncertainty associated with random parameters, one can employ various uncertainty quantification tools. Monte Carlo Simulations (MCS) is a popular approach and has been applied to obtain spatially-distributed probabilistic prediction from the stochastic kinematic wave model [21, 19]. It computes the mean and variances of system states by solving the deterministic governing equations multiple times using a number of realizations of the input random parameters or fields. Statistical moments can be then used as a forecast of the system average response and/or a measure of associated prediction error. However, if one is interested in the distribution tail of system states, a common concern in the risk assessment of many environmental applications [27], the MCS framework becomes expensive due to the requisite large number of simulations for accurate estimation of these higher-order moments.

Another popular uncertainty quantification framework is the generalised polynomial chaos (gPC) expansion method [38] that aims to construct a surrogate relationship between system states and random parameters. Based on the seminal work on Hermite polynomial chaos [10], the gPC method is mathematically robust and is numerically easy to implement with stochastic Galerkin method or various stochastic collocation schemes [16], such as sparse grid method [37], pseudospectral method [39], compressive sensing method and various adaptive schemes [6, 31, 7]. However, the “curse of dimensionality” is a major concern that incurs heavy computational cost when the uncertain random fields have weak two-point correlation functions as happens with, for example, white noise [35, 36].

In recent years the so-called CDF method was proposed to address the parametric uncertainty of kinematic wave models [33]. Based on the concept of (one-point, one-time, Eulerian) fine-grained CDF, it is an extension of early works on the PDF method [20, 23, 28] and has found applications in two-phase flow in porous media [34, 3, 4] and others. Compared to direct simulations of the kinematic wave model, the CDF method offers three major advantages: (1) linearity of new fine-grained CDF equation can drastically reduce the computational cost; (2) one needs only the first statistical moment of solution to the fine-grained CDF equation and (3) its convergence rate is independent of the number of random parameters. Previous studies have since proposed semi-analytical solutions to special cases of the kinematic wave model [33], but a comprehensive and efficient numerical scheme is yet developed for the fine-grained CDF equation. This is the focus of our study.

In this paper, we present a numerical scheme to implement the CDF method for the kinematic wave equation. In Section 2, we provide a brief review of the kinematic wave model and then derive its corresponding fine-grained CDF equation using the CDF method. A comprehensive numerical scheme is then proposed and verified with a simple case in Section 4. We investigate the robustness and salient features of the numerical scheme in Section 4 by comparing the CDF solutions of one kinematic wave system, the Saint-Venant equation, with those obtained from direct MCS of the original kinematic wave model. Key conclusions are drawn in Section 5.

2 Problem Formulation

In this Section we formulate the general problem of stochastic kinematic wave (Section 2.1) and present the CDF formulation from previous work [33] (Section 2.2). Without specifying otherwise, our work is conducted on a complete probability space $(\Omega,\mathcal{F},\mu)$, where $\Omega$ is the collection of events, $\mathcal{F}$ represents a $\sigma$-algebra on sets of $\Omega$, and $\mu$ denotes a probability measure on $\mathcal{F}$. We assume that all random variables have finite second moment. Vectors are represented by lowercase boldface letters, e.g., $\mathbf{q}$, and scalars by lowercase normal letters, .e.g., $k$. The superscript in $\mathbf{q}^{\mathrm{T}}$ represents the vector transpose operator. We reserve $z$ and all its variations ($\boldsymbol{z},Z,\boldsymbol{Z}$) for random variables.

2.1 Kinematic wave equation

Kinematic wave models are popular tools in the study of flow motions driven by gravity and pressure. These models are constructed by postulating functional relationships between the quantity of interest $k(\mathbf{x},t)$ (quantity per unit distance) and its flux $\mathbf{q}(\mathbf{x},t)$ (quantity passing a given point in unit time and distance),

[TABLE]

The kinematic wave model combines the conservation laws of mass and momentum and describes the wave phenomenon with the continuity equation alone:

[TABLE]

subject to the following initial and boundary conditions:

[TABLE]

Here, $S(\mathbf{x},t)$ is a source/sink term, $\mathbf{z}\equiv(z_{1},\dots,z_{\mathrm{N}})$ are $N$ random parameters in the constitutive relationship (1a), and $\mathcal{D}$ is a spatial domain.

In general, the constitutive relationship (1a) and its parameterization $\mathbf{z}$ are often based on physical interpretations of the underlying process. For example, in overland flows one may use Darcy-Weisbach, Chézy, or Manning formulae to describe laminar, turbulent or transitional flow regime, respectively [25].

In many environmental applications, source/sink, functional parameters $\mathbf{z}$, boundary and initial conditions may be subject to epistemic uncertainty. In the case of open-channel flow, $S(\mathbf{x},t)$ represents the rainfall rate and/or inflow/outflow at the tributaries and/or runoff from the ambient terrain, all of which exhibit heterogeneity at various spatio-temporal scales. Meanwhile, past data analysis ([5, 8, 9, 18] and the references therein) of the two functional parameters, e.g. the channel slope and surface roughness coefficient, suggest site-specific statistical distributions to capture their spatial variability. Although data acquisition continues to improve, ubiquitous data sparsity and measurement/interpretation errors render overland flow predictions inherently uncertain. This predictive uncertainty is routinely mentioned as one of the fundamental challenges in flood forecasting [17].

In subsequent analysis, we model the parametric uncertainty as (un/correlated) random fields. In other words, a quantity of interest $k(\mathbf{x},t,\omega)$ varies not only in the physical domain, $(\mathbf{x},t)\in\mathcal{D}\times[0,\infty)$, but also with respect to the probability event $\omega\in\Omega$.

2.2 CDF method

To solve the stochastic kinematic wave equation (1) and obtain a statistical description of $k$ at any space-time point, one can employ the CDF method [33, 34, 3]. In the deterministic setting, we introduce the concept of (one-point, one-time, Eulerian) fine-grained cumulative density function (CDF) of $k$:

[TABLE]

where $\mathcal{H}(\cdot)$ is the Heaviside step function and $K$ is a deterministic value (outcome) that the random quantity $k$ can take at a space-time point $(\mathbf{x},t)$. The Heaviside function is defined here as

[TABLE]

In general, the relation (1a) with random variables $\boldsymbol{z}$ induces randomness in the solution. In this case the function $\Pi$ is random, a function of $\omega$, but in the sequel we continue to refer to it as a “CDF”, consistent with previous literature. With $f_{\mathrm{k}}(K;\mathbf{x},t)$ denoting the single-point probability density function (PDF) of $k(\mathbf{x},t,\omega)$, at fixed $(\boldsymbol{x},t)$, the expectation of the $\Pi$ as defined in (2) yields the single-point CDF $F_{\mathrm{k}}(K;\mathbf{x},t)$,

[TABLE]

Here we emphasize again, that the “CDF” $\Pi$ defined in (2) is an indicator function, and is random. The actual PDF and CDF of $k$ are $f_{\mathrm{k}}$ and $F_{\mathrm{k}}$, respectively, and are deterministic. Following earlier works [33], we multiply the kinematic wave equation (1b) with $\partial\Pi/\partial K$ and obtain a linear fine-grained CDF equation (see Appendix A):

[TABLE]

whose initial and boundary conditions are derived from the physical space relations (1c)&(1d):

[TABLE]

The CDF formulation offers four major advantages. First, comparing to the nonlinear kinematic wave equation (1b), the linear fine-grained CDF equation (4) is easier to solve, albeit in a higher spatial dimension. Second, by invoking (3), one only needs to compute the ensemble average of $\Pi$ to obtain the full statistical distribution, bypassing the need to compute high-order moments. Third, the CDF formulation does not impose any prior assumption on the number of random parameters or on their correlation structures. Lastly, the fine-grained CDF boundary condition (5c) is formulated in a straightforward and unambiguous way, whereas its counterpart in the PDF methods [20, 28], e.g. the probability density of a specific condition, may be hard to formulate or is generally unknown.

3 Numerical scheme for the CDF method

In this section, we present an efficient numerical scheme (Section 3.1) for the CDF formulation to obtain a complete (single space-time point) probabilistic description of $k(\mathbf{x},t)$ via its cumulative density function $F_{\mathrm{k}}(K;\,\mathbf{x},t)$. We then consider a few examples, namely, one-dimensional flow, three-dimensional flow, coupled system of one-dimensional flow and Burgers’ equation, to examine the accuracy and illustrate some salient properties of the numerical scheme (Section 3.2).

3.1 Numerical scheme

We propose an efficient numerical framework that focuses on the computation of the fine-grained CDF equation (4). Due to its linearity, we can use the well-known method of characteristics to solve the $\Pi$ equation. Specifically, we obtain the family of characteristics,

[TABLE]

along which the fine-grained CDF equation (4) reads:

[TABLE]

The original problem (4) is therefore recast to a set of ordinary differential equations (6) and (7), which can be then solved by standard numerical ODE solvers. We employ a third-order Total Variation Diminishing (TVD) Runge-Kutta scheme [11].

By invoking (3), we can create an ensemble of $M$ realization of $\boldsymbol{z}$ in (1a), compute the associated ensemble of solutions $\Pi$ by solving (4), and finally compute $F_{\mathrm{k},M}(K;\boldsymbol{x},t)$, which is an approximation to integral $F_{\mathrm{k}}(K;\mathbf{x},t)$ from the ensemble of its $\Pi$ solutions.

[TABLE]

where $\boldsymbol{z}_{i}$ represents the $i$-th realization of $\boldsymbol{z}$, $\Pi_{i}$ is the corresponding fine-grained CDF solution, $F_{\boldsymbol{z}}(\boldsymbol{z}_{i}^{\prime})$ is the joint probability of random inputs for that realization while $F_{\boldsymbol{z}}(\boldsymbol{z}_{0}^{\prime})=0$ and $F_{\boldsymbol{z}}(\boldsymbol{z}_{M+1}^{\prime})=1$, respectively. For independent random variables $\boldsymbol{z}$, the joint probability $F_{\boldsymbol{z},i}$ would be a tensor product of each random inputs and various collocation techniques can be implemented. For correlated random inputs, one may employ the latest work [12] to compute the high dimensional integral with designed quadrature. We hope to address this subject in future works.

The overall procedure can be summarized as follows:

Generate $M$ realizations of the random parameters $\boldsymbol{z}$. 2. 2.

Solve the fine-grained CDF equation (3) $M$ times via the method of characteristics, one for each realization. 3. 3.

Compute the ensemble average of the $\Pi$ solutions using (8).

In essence, our numerical framework is based on two key features of the CDF formulation: linearity of the fine-grained CDF equation (4) and the ensemble average of its solution is the CDF of system state (3). To implement the method of characteristics, one can take advantage of the large body of literature on the numerical methods for ordinary differential equations.

Alternatively, one may employ the Reynolds decomposition: $\mathcal{A}=\langle\mathcal{A}\rangle+\mathcal{A}^{\prime}$, to represent the random parameters as the sum of their ensemble means $\langle\mathcal{A}\rangle$ and zero-mean fluctuations about the mean $\mathcal{A}^{\prime}$. Then, by taking the ensemble average of the linear fine-grained equation (4), a deterministic equation ensues:

[TABLE]

where $\mathbf{v}_{\mathrm{eff}}$ and $\mathbf{D}$ are the effective velocity and the eddy-diffusivity tensor, respectively. Similar to standard PDF methods [29, 28, 32, 30], such numerical approach (9) requires a closure approximation, such as the large-eddy-diffusivity closure, and its solution is asymptotically exact when $F_{\mathrm{k}}$ varies slowly with $\mathbf{x}$, $K$ and $t$, relative to $\mathbf{v}$ [13, 2].

3.2 Numerical examples

In this section, we present several numerical examples to illustrate the applicability and efficiency of our method.

3.2.1 One-dimensional flow

Let us consider the following kinematic wave model:

[TABLE]

whose corresponding fine-grained CDF equation is:

[TABLE]

Deterministic case

We first examine the accuracy of our numerical scheme by considering deterministic source term and initial and boundary conditions:

[TABLE]

which yields an explicit closed-form solution:

[TABLE]

Now we implement our numerical framework and compute the fine-grained CDF equation (11) with a third-order Runge-Kutta scheme (Appendix B). We compared the numerical results with the analytical solution (13). The error is measured

[TABLE]

where $k(x_{i},t)$ is the solution computed using the method of characteristics, and $x_{i}$ is an $N$-point equidistant grid on $[0,2]$.

Table 1 shows the comparison results for the solution at time $t=0.1$ and one can see that, with a diminishing time-step (from $0.1$ to $0.0125$), the error (in $l_{2}$-norm) $\epsilon$ between our numerical result and the exact solution diminishes to $1.71\times 10^{-5}$. We observe third-order convergence in time for our approach, which matches expectations for this third-order Runge-Kutta scheme.

Stochastic case

We now consider a random source but maintain deterministic initial and boundary conditions,

[TABLE]

where $z$ is a random variable with lognormal distribution:

[TABLE]

whose mean and variance are: $\mu=0$, $\sigma^{2}=0.1$.

The exact solution of the kinematic wave (10) is:

[TABLE]

whose exact CDF can be found as:

[TABLE]

Using the proposed numerical scheme, we solve the stochastic fine-grained CDF equation (11) for a number of realizations of $z$. For each realization, it is found that the errors between exact and numerical solution are in the same order, e.g. $\sim O(10^{-5})$, as shown in the deterministic case. We then compute the ensemble average $F_{\mathrm{k},M}$ based on  (8) . In Fig. 1 (a), the numerical CDF solution at a single space-time point, $F_{\mathrm{k},M}(K;x=0.2,\,t=1)$, shows a good agreement with the analytical solution (18). For a closer examination in Fig. 1 (b), we see that as the realizations number increases, the error ($\epsilon$) rapidly converges and then saturates after $100$ realizations. The error $\epsilon$ here for $F_{\mathrm{k}}$ is defined as:

[TABLE]

where $K_{i}$ is an $N$-point equidistant grid in $K$.

3.2.2 Three-dimensional flow

We now consider a three-dimensional problem:

[TABLE]

where $z$ is a random variable with lognormal distribution: $\ln z\sim\mathcal{N}(0,0.1^{2})$. One can find an exact solution to the three-dimensional flow as follows:

[TABLE]

In Fig. 2, we plot the $l_{2}$-norm error (19) of the CDF solutions $F_{\mathrm{k}}(K;x_{1}=1.3,x_{2}=1.3,x_{3}=1.3,t=1)$ between the exact solution (21) and the one obtained from the proposed CDF scheme, and that from the Monte Carlo Simulations (MCS) of the original system (3.2.2), respectively, at different realization number $M$. We found that the CDF method exhibits excellent accuracy with a fast convergence rate.

3.2.3 Coupled system of one-dimensional flow

We now examine the following system of one-dimensional flow:

[TABLE]

subject to the random initial conditions:

[TABLE]

where $z$ is a random variable with lognormal distribution: $\ln z\sim\mathcal{N}(0,0.1^{2})$.

Its exact solution can be found as

[TABLE]

We introduce two fine-grained CDFs:

[TABLE]

where the variables: $v_{1}=(k_{1}+k_{2})/2$ and $v_{2}=(k_{1}-k_{2})/2$ can help decouple the original system. Now one can employ the proposed CDF scheme to obtain those two marginal fine-grained CDFs for each realization and then recover the fine-grained CDFs of original system states:

[TABLE]

With a number of realizations, the marginal cumulative density functions $F_{\mathrm{k}_{1}}$ and $F_{\mathrm{k}_{2}}$ can be obtained with the approximation formula (8). The $l_{2}$-norm errors (19) between the exact solutions (26) and those from the CDF method, and MCS of the original coupled system, respectively, are plotted in Figure. 3 for $F_{\mathrm{k}_{1}}(k_{1}^{\prime};x=0.3,t=1)$ and $F_{\mathrm{k}_{2}}(k_{2}^{\prime};x=0.3,t=1)$ at different realizations number $M$. It is clear that the CDF scheme provides better accuracy for the same number of realizations and achieves faster convergence rate.

3.2.4 Burgers’ equation

Finally, let us consider a nonlinear example, the Burgers’ equation:

[TABLE]

subject to a stochastic initial condition with lognormal random variable $z$: $\ln z\sim\mathcal{N}(0,0.1^{2})$.

[TABLE]

The initial condition of the Burgers’ equation may lead to shocks at later time. To address such issue, we follow earlier works [34, 1] and employ the Rankine-Hugoniot condition [26] to determine the shock location $x_{\mathrm{s}}$ at each realization. Here the fine-grained CDF solution would be divided by two parts: behind $x_{\mathrm{s}}$, the $\Pi$ satisfies the governing equation whereas it remains its initial condition ahead of $x_{\mathrm{s}}$:

[TABLE]

Following the proposed CDF scheme, we calculate the cumulative density function $F_{\mathrm{k}}$ from solutions of $\Pi$. Figure 4 presents the $l_{2}$-norm errors (19) between the converged solution $F_{\mathrm{k}}(K;x=0.4,\,t=1)$ and those from CDF method, and the MCS simulations of the Burger’s equation, respectively, at different realization numbers $M$. Here the reference solution is obtained from $2000$ MCS simulations.

4 Results and discussion

In this section, we investigate the robustness of our numerical scheme through an environmental application of the kinematic wave model. To be specific, we consider a one-dimensional Manning open-channel flow under three source cases, namely: no source ($S=0$), constant source ($S=\mbox{const}$) and spatially-dependent source ($S=S(x)$). Solutions from our CDF approach and those from the direct simulations would be compared and analyzed.

The Saint-Venant equation is the unidirectional form of shallow water equations and often used to describe flood wave propagation. Using the Manning constitutive relationship, it can be written as:

[TABLE]

where $k(x,t)$ $[\mathrm{m}^{2}]$ is the cross-sectional area of a channel occupied by the fluid at a point $x$ along the channel length, $q(x,t)$ $[{\mathrm{m}}^{3}/{\mathrm{s}}]$ describes the volumetric flow rate, $S$ $[\mathrm{m}^{2}/\mathrm{s}]$ is the lateral inflow rate of tributaries and/or upstream rainfall rate, $s_{0}(x)$ is the slope of the channel bed, and $C_{\mathrm{M}}(x)$ [s/m1/3] denotes the Manning’s roughness coefficient. We consider $C_{\mathrm{M}}$ and $s_{0}$ as random fields. Thus, we have the flux function formulation (1a), with $\boldsymbol{z}=(z_{1},z_{2})=(C_{\mathrm{M}},s_{0})$.

The Saint-Venant model (37) provides a good approximation of flood waves for a Froude number smaller than one, in which the main disturbance is carried downstream only by the kinematic waves while dynamic waves (long gravity waves) attenuate rapidly [14]. We also note that the kinematic wave model neglects influence on the river upstream of the junction [14] and the backwater effects (upstream propagation caused by local acceleration, convective acceleration and pressure), the flow rate throughout the flow domain is therefore non-negative, $q(x,t)\geq 0$. Without loss of generality, the approach presented below can also incorporate higher physical dimensions and other types of constitutive relationship, such as the Chézy formula, to represent a balance between the friction at the channel bottom and the gravitational force.

With no analytical solution available, the Saint-Venant equation (37) is often solved numerically in the form of flux $q$ [33]:

[TABLE]

Such equation would be directly solved via a fifth-order weighted essentially non-oscillitorary (WENO) scheme in space and third-order TVD Runge-Kutta scheme (Appendix C).

In order to compare with the direct simulation of the hyperbolic equation above (37), we apply the CDF formulation to solve for flux $q$. The corresponding fine-grained CDF equation of flux, $\Pi=\mathcal{H}\left[Q-q(x,t)\right]$ is:

[TABLE]

Following statistical data from earlier analysis [5, 8, 9, 18], both slope $s_{0}(x)$ and the Manning coefficient $C_{\mathrm{M}}(x)$ are treated as stationary random fields exhibiting lognormal distribution and exponential correlation function with correlation length $\lambda$. The mean and standard deviation are $0.01$, $0.0025$ and $0.037$ [s/m1/3], $0.00925$ [s/m1/3], for $s_{0}(x)$ and $C_{\mathrm{M}}(x)$, respectively. We employ Gaussian process model to produce $M=1000$ realizations for each random parameters. Flow rate at initial time and the inlet is set as $q_{\mathrm{in}}(x,t=0)=0.5$ [m3/s] and $q_{\mathrm{bx}}(x=0,t)=\mathrm{max}\left(\sin\pi t,0.5\right)$ [m3/s], respectively.

We conduct $1000$ MCS simulations for the original kinematic wave equation (37) and use it as the reference solutions. Figure 5 illustrates the cumulative density function of flux $F_{\mathrm{q}}(Q;x=1,\,t=1)$, under three source conditions: $S=0$, $S=1$ and $S=x$. It is clear that solutions from the CDF method match well with those from the direct MCS simulations. For such desired ensemble accuracy, the fine-grained CDF equation at each realization requires significantly less CPU time to compute (Table 2). This is one of the main advantages of our approach using characteristics: although using direct solvers for the kinematic wave equation (in this case WENO-Roe finite-difference schemes) seems straightforward, using the method of characteristics by first deriving the $\Pi$ equation results in more than 50X speedup with a negligible compromise in accuracy.

We now examine the convergence rates of the CDF method and direct MCS. Their relative error (19) to the benchmark MCS solution $F_{\mathrm{q}}(Q;x=1,\,t=1)$ are plotted in Fig. 6 for different realization numbers $M$. We find that the proposed CDF scheme yields smaller errors than those from the MCS for the same number of realizations and overall exhibits faster convergence rate.

5 Conclusion

In this paper, we presented a numerical framework to implement the CDF method in order to obtain the cumulative density function (CDF) of the system state in the kinematic wave model. The approach relies on solving the fine-grained CDF equation derived from the CDF method. Its accuracy and robustness were investigated via comparison with direct MCS of several numerical examples and a kinematic wave system, the one-dimensional Saint-Venant equation. Our analysis leads to the following conclusions:

• •

In contrast to previous CDF approach that focuses on the derivation and computation of an ensemble-averaged equation of system state distribution, we directly solve the fine-grained CDF equation, which is exact from the original stochastic system.

• •

At each single realization, the proposed numerical scheme proves to be computationally more efficient in solving the linear fine-grained CDF equation than the direct simulation of the nonlinear kinematic wave equation.

• •

In obtaining the cumulative density function of the system state, our scheme exhibits superior convergence rate and requires fewer realizations than the direct MCS.

• •

Faster convergence rate can be achieved with improved quadrature rules and would be the focus of future works.

Appendix A Derivation of Raw CDF Equation

Following the definition of $\Pi$ in (2), we find its spatial and temporal derivatives as

[TABLE]

For smooth solutions, we multiply the kinematic wave equation (1b) with $\partial\Pi/\partial K$ and substitute the derivates above (39):

[TABLE]

where the operator $\nabla_{\mathbf{z}}(\cdot)$ is defined as:

[TABLE]

We note here that derivative of the fine-grained CDF in the probability space is the dirac delta function, e.g. $\partial\Pi/\partial K=\delta(K-k)$. Using its sifting property, $g(k)\delta(K-k)=g(K)\delta(K-k)$ for any test function $g(\cdot)$, all $k(\mathbf{x},t)$ in the equation (40) can be replaced by $K$ and leads to a linear fine-grained CDF equation (4). Interested readers can refer to earlier studies [20, 33] for detailed derivations.

Appendix B TVD Runge-Kutta scheme for characteristic lines

The Runge-Kutta method is a simple and robust iterative scheme. Here we employ a third-order TVP Runge-Kutta scheme [11] to solve the ordinary differential equations (6) and solution at the $h$-th stencil and the $n$-th time-step, $k_{h}^{n}$, can be obtained in a three-step iteration:

[TABLE]

Appendix C WENO-Roe Numerical scheme for kinematic wave equation

The WENO scheme [22] provides a nonlinear adaptive procedure to automatically select the smoothest local stencil in numerical approximation of fluxes. We first take the change of variable, $q^{\prime}=(q\,C_{\mathrm{M}}/\sqrt{s_{0}})^{3/4}$, and rewrite the one-dimensional kinematic wave equation (37) as:

[TABLE]

The spatial domain $x\in[a,b]$ is divided into $N_{\mathrm{x}}$ stencils of size $\Delta x$ and let us denote $q^{\prime}_{h}$ as the value of $q^{\prime}$ at the $h$-th stencil. Using a fifth-order finite difference WENO-Roe scheme, we can numerically approximate the partial differential equation with an ordinary differential equation:

[TABLE]

and $u(q^{\prime})=\sqrt{s_{0}}q^{\prime 4/3}/C_{\mathrm{M}}$. The superscripts ± refers to the direction from which one interpolates the half stencils, in other words, if $a_{h+\frac{1}{2}}>0$, the value $u(q^{\prime})^{-}_{h-\frac{1}{2}}$ is approximated from the upwind direction. Exact expressions of those interpolations can be found in earlier studies [22].

Now we can solve the ordinary differential equation with a third-order TVD Runge-Kutta scheme (Appendix B). Boundary conditions at the $h_{-2}$-th and $h_{-1}$-th stencils are taken as the same value as that at $h_{0}$-th.

We validated our WENO-Roe scheme via the deterministic example in section 3.2 with a time step $\Delta t=0.0001$. Table 3 shows the comparison between our numerical result and the analytical solution (13) at time $t=0.1$: as the space step $\Delta x$ drops (from $0.05$ to $0.00625$), the 2-norm difference ($\epsilon$) between the two solutions is reduced to $1.97\times 10^{-9}$ while the WENO-Roe convergence rate is approximately $4.5$, providing reasonable validation of the fifth-order approach.

Acknowledgment

P. Wang and M. Cheng were partially funded by the National Natural Science Foundation of China (Grant No. 11571028), the National Key Research and Development Program of China (Grant No. 2017YFB0701702) and the Recruitment Program of Global Experts. A. Narayan was partially supported by AFOSR FA9550-15-1-0467. X. Zhu was partially supported by Simons Foundation.