Uncertainty quantification of floodplain friction in hydrodynamic models

This study proposes a framework to estimate the uncertainty of hydrodynamic models on floodplains. The traditional floodplain resistance formula of Pasche (1984) (based on Lindner, 1982) used for river modelling as well as the approaches of Baptist et al. (2007), Järvelä (2004), and Battiato and Rubol (2014) was considered for carrying out an uncertainty quantification (UQ). The analysis was performed by means of three different methods: traditional Monte Carlo (MC), first-order second-moment (FOSM) and metamodelling. Using a two-dimensional hydrodynamic model, a 10 km reach of the River Rhine was simulated. The model was calibrated with water level measurements under steady flow conditions and then the analysis was carried out based on flow velocity results. The compared floodplain friction formulae produced qualitatively similar results, in which uncertainties in flow velocity were most significant on the floodplains. Among the tested resistance formulae the approach from Järvelä (2004) presented on average the smallest prediction intervals, i.e. the smallest variance. It is important to keep in mind that UQ results are not only dependent on the defined input parameter deviations, but also on the number of parameters considered in the analysis. In that sense, the approach from Battiato and Rubol (2014) is still attractive for it reduces the current analysis to a single parameter, the canopy permeability. The three UQ methods compared gave similar results, which means that FOSM is the less expensive one. Nevertheless it should be used with caution as it is a first-order method (linear approximation). In studies involving dominant non-linear processes, one is advised to carry out further comparisons.

1 Introduction

Flow resistance can be considered as the contribution of four components, according to Rouse (1965): (a) surface, (b) form, (c) wave and (d) flow unsteadiness resistances. Not only that, it is a complex phenomenon dependent on Reynolds number, relative roughness, cross-sectional geometry, channel non-uniformity, Froude number and degree of flow unsteadiness. Also, Yen (2002) affirms that flow resistance interacts “in a non-linear manner such that any linear separation and combination is artificial”.

There are several approaches available in the literature for determining flow resistance coefficients of vegetated floodplains in numerical models. These approaches are basically divided into four categories: rigid or flexible, and emergent or submerged vegetation. They aim to determine the resistance exerted by the vegetation on the flow based on physical properties such as vegetation height and width, stem diameter and density, etc. Recent research on flow resistance of emergent floodplain vegetation is given in Aberle and Järvelä (2013) and a review of vegetated flow models can be found in Nikora et al. (2007). For an overview of the main vegetation friction laws available the reader is referred to the review given in Shields et al. (2017).

Even though much work has been done in applying different approaches to include vegetation-induced resistance effects in hydrodynamic calculations, the majority of these studies were verified only under laboratory-scale conditions. A gap between those results and river engineering projects still exists. While free surface information on flooded areas can be well approximated from river channel measurements, flow velocity cannot. And because floodplain measurements are usually not available, model performance is neglected at those areas. That means, when flood scenarios belong to the scope of a project or study, attention should be given to this matter. A way to address this problem is to consider a probabilistic approach and to carry out an uncertainty quantification (UQ) of the floodplain friction. Uncertainty in the context of fluid dynamics is defined as a potential deficiency of the simulation process, according to Walters and Huyse (2002). Straatsma and Huthoff (2010) considered floodplain friction parametrization to be an important source of uncertainty. Also, Di Baldassarre et al. (2010) compared and discussed deterministic and probabilistic approaches for floodplain mapping. They concluded that due to uncertainties related to flood-event statistics the probabilistic approach was considered to be a more correct representation.

Some studies can be found in the literature that involve the estimation of uncertainty related to floodplains and the resistance coefficient. Apel et al. (2004) presented a flood risk assessment by means of a simple hydrological flood routing model in the Lower Rhine applying a Monte Carlo (MC) framework. Pappenberger et al. (2005) conducted an analysis using a one-dimensional hydraulic model with a generalized likelihood uncertainty estimation. Their results showed that many parameter sets (channel and floodplain) can perform equally well even with extreme values. Brandimarte and Di Baldassare (2012) showed that the deterministic approach underestimates the design flood profile in hydraulic modelling and proposed an alternative approach based on the use of uncertain flood profiles. Altarejos-García et al. (2012) used the point-estimate method for carrying out uncertainty quantification as an alternative to MC approaches to get estimates of the mean and variance of water depth and velocity. They considered the roughness coefficient as the main source when assessing the uncertainty in river flood modelling. Domeneghetti et al. (2013) proposed a methodology to derive probabilistic flood maps, taking into account several sources of uncertainty. Willis et al. (2016) concluded that hydrodynamic modelling can be improved by increasing the number of frictional surfaces; however, they draw attention to the numerical scheme choice, which might lead to much larger errors.

In this context, a framework to estimate the uncertainty of hydrodynamic models on floodplains due to vegetation is proposed in the current study. Within a large number of different floodplain friction methods there is still no generally accepted choice for large-scale applications. Thus, the outcomes of the uncertainty quantification can assist to identify a well-suited friction method for practical use. A two-dimensional hydrodynamic model is calibrated with floodplain friction formulations, with which input parameter uncertainties are associated based on a practical range of variation. After defining the variation ranges for sensitive input parameters, the UQ is carried out with different methods for comparison. In other words, the quantification of uncertainties in flow velocity simulation is addressed by considering uncertainties in floodplain friction parameters. In the next section four chosen floodplain resistance formulae are described and analysed. Then the concept of uncertainty quantification is briefly explained and three different methods are presented in the third part. The fourth section provides information on the case study including a brief description of the hydrodynamic model, parameters used for model calibration and the defined uncertainties for the selected input parameters needed for carrying out the analysis. In the fifth section results are presented and discussed, from which conclusions are drawn in the last part of the paper.

2 Floodplain friction

Vegetation found on river banks and floodplains plays an important role in the flow velocity profile and, therefore, on hydraulic roughness. Current research aims to relate vegetated floodplain properties to their “hydraulic signatures” and to incorporate the complex nature of vegetation characteristics into floodplain friction models. According to Shields et al. (2017), there are no established practices for defining flow-dependent vegetation roughness values and incorporating them into hydrodynamic models. Additionally, model calibration is usually carried out with measurements taken in the main channel, and seldom (if ever) on floodplains. Thus, model response on floodplains cannot be verified and only relative conclusions can be made. It is under these circumstances that UQ is especially useful for quantifying the probability of results. Basically the available approaches for vegetation friction formulation are subdivided into emergent or submerged and rigid or flexible.

For the current study four out of several vegetation friction formulations (see Aberle and Järvelä, 2015; Shields et al., 2017) are considered for floodplains: Lindner (1982) and Pasche (1984), Baptist et al. (2007), Järvelä (2004), and Battiato and Rubol (2014). The first approach is a recommended practice by the German Association for Water, Wastewater and Waste (DVWK, 1991) for hydraulic calculations and it is commonly used in the German Federal Waterways Engineering and Research Institute's (BAW) projects. The second and third approaches represent the rigid (Baptist et al., 2007) and flexible (Järvelä, 2004) approximations. Lastly, the approach from Battiato and Rubol (2014) is chosen as it proposes a completely different concept based on porous medium flow.

2.1 Lindner and Pasche

The modified formulation from Pasche (1984), based on Lindner (1982), was developed for rigid emergent vegetation. The Darcy–Weisbach friction factor for vegetation (f_v) can be obtained after the bulk drag coefficient (C_D) is iteratively calculated by the following equations:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& U=\sqrt{{\displaystyle \frac{\mathrm{8}g{R}_{\mathrm{h}}{S}_{\mathrm{0}}}{f}}},\mathrm{with}f=f_{\mathrm{b}}+f_{\mathrm{v}},\\ & {\left({\displaystyle \frac{U_i}{U}}\right)}^{\mathrm{2}}=\mathrm{1.15}{\left(W_{\mathrm{l}}\sqrt{m}\right)}^{-\mathrm{0.48}}+\mathrm{0.5}{\left(W_{\mathrm{w}}\sqrt{m}\right)}^{\mathrm{1.1}},\\ & \mathrm{0.03}=\mathrm{0.9}{\left({\displaystyle \frac{W_{\mathrm{l}}}{C_{\mathrm{D}\mathrm{1}}}}\right)}^{-\mathrm{0.7}}{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}g{W}_{\mathrm{l}}{S}_{\mathrm{0}}}{U^{\mathrm{2}}}}\right)}^{-\mathrm{1.5}},\\ & W_{\mathrm{w}}=\mathrm{0.24}{W}_{\mathrm{l}}^{\mathrm{0.59}}(C_{\mathrm{D}\mathrm{1}}D{)}^{\mathrm{0.41}},\\ & \mathrm{Fr}={\displaystyle \frac{U}{\sqrt{gh}}}\simeq {\displaystyle \frac{y^{*}(y^{*\mathrm{2}}-\mathrm{1})}{\mathrm{2}[y^{*}-\mathrm{1}/(\mathrm{1}-D\sqrt{m}\left)\right]}},\\ & C_{\mathrm{D}}=\mathrm{1.31}{\left({\displaystyle \frac{U_i}{U}}\right)}^{\mathrm{2}}+\mathrm{2}{\mathrm{Fr}}^{-\mathrm{2}}(\mathrm{1}-y^{*}),\\ & f_{\mathrm{v}}=\mathrm{4}{C}_{\mathrm{D}}mDh,\end{array}\end{array}$$

where R_h is the hydraulic radius, S₀ is the bottom slope, f_b is the bottom friction, U is the approach velocity (upstream), U_i is the calculated velocity (downstream), W_l and W_w are the wake length and width respectively, C_D1 is the drag coefficient for a single stem, m is the number of stems per square metre, D is the stem diameter, h is the water depth, and g is the gravitational acceleration.

2.2 Baptist et al.

The approach from Baptist et al. (2007) was developed for rigid vegetation. They modelled the vegetation resistance force as the drag force on an array (random or staggered) of rigid cylinders with uniform properties. The velocity profile is calculated for two conditions: non-submerged (emergent) and submerged vegetation. For the case of emergent vegetation a uniform velocity is assumed. For the case of submerged vegetation the velocity profile is subdivided into a uniform velocity zone (within the vegetation) and logarithmic velocity zone (above the vegetation). Both conditions combined, after some algebra and use of genetic programming, give the following expression for the Chézy coefficient induced by bottom and vegetation friction (C):

$$C={\left({\displaystyle \frac{\mathrm{1}}{C_{\mathrm{b}}^{\mathrm{2}}}}+{\displaystyle \frac{C_{\mathrm{D}}mDh}{\mathrm{2}g}}\right)}^{-\mathrm{0.5}}+{\displaystyle \frac{\sqrt{g}}{\mathit{\kappa }}}\ln \left({\displaystyle \frac{h}{H}}\right),$$

where H is the vegetation height, C_b is the Chézy coefficient of the bed and κ is the von Kármán constant. The corresponding Darcy–Weisbach friction factor can be then obtained by

$$\begin{array}{}\text{(2)}& f={\displaystyle \frac{\mathrm{8}g}{C^{\mathrm{2}}}}.\end{array}$$

2.3 Järvelä

The approach from Järvelä (2004) was developed for flexible vegetation. It is based on the leaf area index (LAI), a dimensionless quantity that characterizes plant canopies. The LAI is defined as the one-sided leaf area per unit projected area in canopies. The Darcy–Weisbach friction factor for vegetation (f_v) can be calculated by the following relation:

$$\begin{array}{}\text{(3)}& f_{\mathrm{v}}=\mathrm{4}{C}_{D\mathit{\chi }}\mathrm{LAI}{\left({\displaystyle \frac{U}{U_{\mathit{\chi }}}}\right)}^{\mathit{\chi }}{\displaystyle \frac{h}{H}},\end{array}$$

where χ is the species-specific vegetation parameter (Vogel exponent), C_Dχ is the species-specific drag coefficient, U is the mean flow velocity, and U_χ is a normalizing value and is defined as the lowest flow velocity used in determining χ. U_χ is usually 0.1 m s⁻¹ and it will be considered constant.

2.4 Battiato and Rubol

The approach from Battiato and Rubol (2014), developed for submerged vegetation, follows the concept of coupling an incompressible fluid flow with a porous medium flow. Although it is conceptually suited for rigid vegetation, this approach was also successfully validated with flexible vegetation (see Rubol et al., 2018). The main advantage of this approach lies in the representation of the drag force by a single parameter, i.e. the canopy permeability (K). The volumetric discharge per unit width through a vegetated channel (Q_w) can be determined from direct integration of the velocity over depth, obtained from the solution of the coupled log–law and Darcy–Brinkman equations:

$$\begin{array}{rl}& L=h-H,\mathit{\delta }={\mathrm{LH}}^{-\mathrm{1}},\mathit{\lambda }=H{K}^{-\mathrm{0.5}},q=\mathit{\rho }g{S}_{\mathrm{0}}{H}^{\mathrm{2}}{\mathit{\mu }}_{\mathrm{t}}^{-\mathrm{1}},\\ & {\mathit{\mu }}_{\mathrm{t}}=\mathit{\rho }{\mathit{\kappa }}^{\prime }H{u}_{*},u_{*}=(g{S}_{\mathrm{0}}L{)}^{\mathrm{0.5}},\\ & C=\mathrm{0.5}\mathit{\delta }{\mathit{\lambda }}^{-\mathrm{1}}\mathrm{csch}\left(\mathit{\lambda }\right),U^{*}={\mathit{\lambda }}^{-\mathrm{2}}+\mathit{\delta }{\mathit{\lambda }}^{-\mathrm{1}}\coth \left(\mathit{\lambda }\right),\\ & Q_{\mathrm{w}}=qH\mathit{\{}{\mathit{\lambda }}^{-\mathrm{2}}+C{\mathit{\lambda }}^{-\mathrm{1}}(e^{\mathit{\lambda }}-e^{-\mathit{\lambda }})+\mathit{\delta }[(\mathrm{1}+\mathit{\delta })\ln (\mathrm{1}+\mathit{\delta }),\\ & +U^{*}-\mathit{\delta }\left]\mathit{\right\}}\end{array}$$

where ρ is the density of water, μ_t is the turbulent viscosity, κ^′ is the reduced von Kármán constant for vegetated channels (${\mathit{\kappa }}^{\prime }=\mathrm{0.19}$) and u_* is the friction velocity. Under emergent conditions (h<H) an approximation is made, in which the velocity profile is considered constant within the canopy. Thus, Q_w is linearized with regard to h by applying $Q_{\mathrm{w}}=Q_{\mathrm{w}}h{H}^{-\mathrm{1}}$. The Darcy–Weisbach friction factor can be then calculated by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathit{\tau }}_H=\mathit{\rho }g{S}_{\mathrm{0}}L,U=Q_{\mathrm{w}}{h}^{-\mathrm{1}},\\ & f={\displaystyle \frac{\mathrm{8}{\mathit{\tau }}_H}{\mathit{\rho }{U}^{\mathrm{2}}}}.\end{array}\end{array}$$

2.5 Overall comparison

From now on the presented floodplain friction formulations will be referred to as LIND, BAPT, JAER and BATT, respectively. The formulae will be analysed in terms of the total Darcy–Weisbach friction factor calculated as $f=f_{\mathrm{b}}+f_{\mathrm{v}}$, with f_b and f_v being the bottom and vegetation friction, respectively. The four expressions are then given by

$$\begin{array}{}\text{(LIND)}& {\displaystyle }f=f_{\mathrm{b}}+\mathrm{4}{C}_{D}mDh,\text{(BAPT)}& {\displaystyle }f=\mathrm{8}g{\left[{\left({\displaystyle \frac{\mathrm{1}}{C_{\mathrm{b}}^{\mathrm{2}}}}+{\displaystyle \frac{C_{\mathrm{D}}mDh}{\mathrm{2}g}}\right)}^{-\mathrm{0.5}}+{\displaystyle \frac{\sqrt{g}}{\mathit{\kappa }}}\ln \left({\displaystyle \frac{h}{H}}\right)\right]}^{-\mathrm{2}},\text{(JAER)}& {\displaystyle }f=f_{\mathrm{b}}+\mathrm{4}{C}_{\mathrm{D}}\mathrm{LAI}{\left({\displaystyle \frac{U}{U_{\mathit{\chi }}}}\right)}^{\mathit{\chi }}{\displaystyle \frac{h}{H}},\text{(BATT)}& {\displaystyle }f=\mathrm{8}g{\displaystyle \frac{S_{\mathrm{0}}H}{U^{\mathrm{2}}}}\left({\displaystyle \frac{h}{H}}-\mathrm{1}\right).\end{array}$$

In LIND and BAPT there is a direct dependency between the term mD and the friction factor. The same analysis is valid for C_D in the first three formulae. The relation h∕H is found in some form in all the approaches which include submerged vegetation. Furthermore, a similar relation between the bottom friction C_b and the friction factor in BAPT is also observed in BATT. While the first three approaches present an explicit term for the bottom friction, in BATT the expression can be rearranged so that a Chézy-like term is found as a function of H.

3 Uncertainty quantification (UQ)

Numerical models represent only an approximation of the observed process. The measured difference between the model and the observation can be considered either as error or uncertainty. Walters and Huyse (2002) defined these two concepts as follows:

• Error: a noticeable lack in the modelling process, not due to a lack of knowledge; (Deterministic)

• Uncertainty: a potential shortcoming in the modelling process due to a lack of knowledge. (Stochastic)

Uncertainty quantification aims to describe the system reliability by combining the uncertainties in the basic components (variables) of the system. The framework of the numerical model used to represent the system characterizes the interactions of the basic components. The overall response of the system is described by the performance function Y:

$$Y=f(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_n)=f\left(\mathit{x}\right),$$

where x is the vector of input variables of the system and n is the number of variables.

The analysis yields the combined effect of all input variables that significantly contribute to the performance function. The results can be represented in terms of reliability or risk. Reliability refers to a prediction interval (PI) of the form [Y_a,Y_b] associated with a known probability $P(Y_a\le Y\le Y_b)$. PI can be expressed as the difference $|Y_a-Y_b|$ with a probability P. Risk refers to the probability of failure (P_f) with respect to a threshold value Y_c, i.e. the probability that Y>Y_c. P_f is directly expressed as the calculated probability.

In the current study Y represents hydrodynamic model results (deterministic) and x contains the input parameters related to each friction formulation. By defining each x_i through a probability distribution based on a potential uncertainty, it is possible to estimate the probability of occurrence of Y. Finally model results are not evaluated as an absolute value any more (e.g. like simulated water depth or flow velocity), but either as a measure given by the absolute difference between Y_a and Y_b with probability P, or as a direct probability P of exceedance of a threshold Y_c.

Three probabilistic methods are chosen for the UQ: first-order second-moment, Monte Carlo (MC) and metamodelling. The first method is based on the method of moments and requires the calculation of the model sensitivities (first-order derivatives). The MC method requires the simulation of a large number of random experiments and is the most expensive in terms of computing time. The metamodelling method is based on random experiments (MC) with the benefit that it requires far fewer samples. Polynomial chaos (PC) is a type of metamodelling technique, which is chosen for the present study. Further details on each method will be given in the following sections.

3.1 First-order second-moment

Moment method approximations are obtained from the truncated Taylor series expansion about the expected value of the input parameter. The first-order second-moment (FOSM) method uses the first-order terms of the series and requires up to the second moments of the uncertain input variables for estimating the output variance of a system. The variance of the performance function σ²(Y) is given by the following:

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}^{\mathrm{2}}\left(Y\right)=\sum _{i=\mathrm{1}}^n\left[{\left({\displaystyle \frac{\partial f}{\partial x_i}}\right)}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}\left(x_i\right)\right].\end{array}$$

It should be noted that the FOSM method is suited as long as (a) the input variables are statistically independent and (b) the linearity assumption is valid, i.e. the first-order approximation is enough to describe the sensitivity of the system. If Y is non-linear, e.g. hydro- and morphodynamic models, one should make sure that the value of σ(x_i) is small. Otherwise, Y might be over- or underestimated. The reader is referred to Dettinger and Wilson (1981), Yen et al. (1986) and Sitar et al. (1987) for further details on FOSM.

3.2 Monte Carlo

Monte Carlo simulation is a probabilistic method in which a very large number of similar random experiments form the basis. An attempt is made to solve analytically unsolvable or complicated solvable problems with the help of probability theory. The law of large numbers makes up one of the main aspects of the method. The random experiments can be carried out in computer calculations in which (pseudo-)random numbers are generated with suitable algorithms to simulate random events.

The basic steps of a MC method can be described as follows:

1. Sample the input random variables x from their known or assumed probability density function N times.

2. Calculate the deterministic output Y for each input sample.

3. Determine the statistics of the distribution of Y (e.g. mean, variance).

Step (2) should be repeated N times, which presents this method's main drawback. Also, the input variables are considered to be statistically independent, otherwise the joint probability distribution is required. The advantage is its robustness, because independently from the nature of Y (linear or non-linear), the method will always deliver reliable results as long as the number of samples (N) is sufficiently large.

3.3 Metamodelling

Metamodelling attempts to offset the increased cost of probabilistic modelling by replacing the expensive evaluation of model calculations with a cost-effective evaluation of surrogates. Polynomial chaos is a powerful metamodelling technique that aims to provide a functional approximation of a computational model through its spectral representation of uncertainty based on polynomial functions. A more detailed introduction to the PC method can be found in Marelli and Sudret (2017).

Spectral-based methods allow for an efficient stochastic reduced basis representation of uncertain parameters in numerical modelling. By means of a truncated expansion to discretize the input random quantities it is possible to reduce the order of complexity of the system. Let us consider the uncertain variable A as a function of the input variable vector x, time t and the random variable vector ξ. A can represent velocity, density or pressure in a stochastic fluid dynamics problem and is defined as a polynomial expansion:

$$A(\mathit{x},t,\mathit{\xi })\approx \sum _{j=\mathrm{0}}^M{a}_j(\mathit{x},t){\mathrm{\Psi }}_j\left(\mathit{\xi }\right),$$

where a_j(x,t) is the deterministic component, Ψ_j(ξ) is the random basis function corresponding to the jth mode and ξ characterizes the uncertainty in the parameter. The polynomial chaos expansion is approximated by a discrete sum taken over the number of output modes M defined as follows:

$$M+\mathrm{1}={\displaystyle \frac{(n+p)\mathrm{!}}{n\mathrm{!}p\mathrm{!}}},$$

where p is the degree of the polynomial and n is the number of random dimensions. The statistics of the distribution for the model output at a specific position and time can be calculated using the coefficients and the basis functions. The mean and variance of the solution is given respectively by

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& E\left[A\right(\mathit{x},t,\mathit{\xi }\left)\right]=\underset{R}{\int }A(\mathit{x},t,\mathit{\xi })w\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi }=a_{\mathrm{0}}(\mathit{x},t),\\ & \mathrm{Var}\left[A\right(\mathit{x},t,\mathit{\xi }\left)\right]=\underset{R}{\int }{\left[A(\mathit{x},t,\mathit{\xi })-a_{\mathrm{0}}(\mathit{x},t)\right]}^{\mathrm{2}}w\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi },\end{array}\end{array}$$

with w(ξ) being the weight function of the polynomial and R its support range. When the input uncertainty is Gaussian (normal) the basis function Ψ(ξ) takes the form of a multi-dimensional Hermite polynomial, so that $R=(-\mathrm{\infty },+\mathrm{\infty })$.

In this study, the non-intrusive polynomial chaos (NIPC) method will be considered. The main objective of this method is to obtain the polynomial coefficients without modifying the original model. This approach considers the deterministic model as a “black box” and approximates the polynomial coefficients based on model evaluations. The advantage is that this method requires much fewer evaluations of the original model (with regard to MC) for providing reliable results (at least 1 order of magnitude). The main disadvantage is that it is an additional approximation in the modelling framework, thus leading to further loss of information of the physical process. The reader is referred to Hosder and Walters (2010) for further details on the application of the NIPC method. The implementation of the method was done in Python with the help of the OpenTURNS package (Baudin et al., 2017).

4 Case study

The current study focuses on a reach of the river Rhine used for numerical tests by the German Federal Waterways Engineering and Research Institute (BAW). It is an 11 km long section of the lower Rhine located between kilometres 738 and 750, near Düsseldorf (Germany). The model topography is based on the digital elevation model provided by the German Federal Centre for Information Technology (2017). The model has been extensively tested and calibrated for a wide spectrum of discharges. A constant discharge of 7870 m³ s⁻¹ was imposed at the upstream boundary and the corresponding free surface at the downstream boundary. These conditions represent a flood scenario with a probability of occurrence larger than HQ₅ (LUA, 2002), i.e. a discharge with an annual probability of occurrence of 1∕5 or likely to be observed within five years (return period of five years). In recent years flood studies are receiving more and more attention as part of BAW's activities. For that reason the current motivation is to understand how sensitive numerical models are to floodplain friction under flood conditions and how this might affect the hydrodynamics of navigation channels. An overview of the study area is presented in Fig. 1, where the red polygon delimits the boundaries of the numerical model.

Figure 1 Lower Rhine river topography and numerical model boundaries (red polygon) nearby Düsseldorf (a), and friction zones definition in the numerical model (b). The flow direction is north. (Satellite image a is the intellectual property of Esri and is used herein under license. Copyright © 2019 Esri and its licensors.)

4.1 Hydrodynamic model

A numerical model is used to simulate the flood scenario. In the BAW, studies carried out in large-scale river projects (10¹–10² km) usually make use of the hydrodynamic model TELEMAC-2D (Galland et al., 1991; Hervouet and Ata, 2017). It is a two-dimensional finite-element software (finite-volume also available) for solving the shallow water equations, a set of partial differential equations derived from the integration of the Navier–Stokes equations over the vertical axis. Thus, the equations for the conservation of mass and momentum in two dimensions should be solved.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial uh}{\partial x}}+{\displaystyle \frac{\partial vh}{\partial y}}=\mathrm{0},\\ & {\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}=-g{\displaystyle \frac{\partial (z_b+h)}{\partial x}}+\mathit{\nu }\left({\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial y^{\mathrm{2}}}}\right)\\ & -{\displaystyle \frac{{\mathit{\tau }}_x}{\mathit{\rho }h}}+S_x,\\ & {\displaystyle \frac{\partial v}{\partial t}}+u{\displaystyle \frac{\partial v}{\partial x}}+v{\displaystyle \frac{\partial v}{\partial y}}=-g{\displaystyle \frac{\partial (z_b+h)}{\partial y}}+\mathit{\nu }\left({\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial y^{\mathrm{2}}}}\right)\\ & -{\displaystyle \frac{{\mathit{\tau }}_y}{\mathit{\rho }h}}+S_y,\end{array}\end{array}$$

where x and y are the Cartesian coordinates, h is the water depth, z_b is the bottom elevation, u and v are the components of the velocity field, ν is the fluid viscosity (which may be constant or given by a turbulence model), τ_x,τ_y are the shear stress components, and S_x and S_y are any additional source term components of momentum (e.g. wind stress, external forces). The bottom shear stress is bound to the depth-averaged velocity by the quadratic law first introduced by Taylor and Shaw (1920):

$$\mathit{\tau }=c_{\mathrm{f}}\mathit{\rho }\mathit{u}\left|\mathit{u}\right|,c_{\mathrm{f}}=c_{\mathrm{fb}}+c_{\mathrm{fv}},$$

The friction coefficient (c_f) is equal to the sum of the bottom friction (c_fb) and the friction due to vegetation (c_fv). The bottom friction can usually be determined by traditional friction laws relating open-channel flow velocity to resistance coefficients (e.g. Manning, Chézy, Darcy–Weisbach, Nikuradse). However, on floodplains the velocity profile strongly depends on the vegetation height and morphology. Thus, specific flow resistance formulae have been developed for determining the vegetation drag (see Sect. 2).

The model consists of an unstructured triangular mesh composed of 56 825 points and 112 360 elements. The resolution varies from about 2.5 m in the main channel to about 30 m on the floodplains and the model mesh covers an area of ca. 8 km². A constant discharge upstream and a constant water level downstream are imposed at the open boundaries, as mentioned above. A time step of 1 s guarantees a Courant number below 1 and it is used to simulate 24 h, which takes about 9 min with the LIND formulation using 160 processors of the BAW's high-performance computing system. The other three formulations for floodplain friction are about three times faster (3.5 min) as there is no iteration step.

This numerical model was extensively investigated from the point of view of sediment transport and morphodynamics (Backhaus et al., 2014; Riesterer et al., 2014), from which the calibrated hydrodynamic model was taken as a starting point for the current study. Water level measurements for a discharge of 7870 m³ s⁻¹ available along the river channel axis were used to calibrate the bottom friction in the numerical model, as a representation of a flood scenario. A root mean square error (RMSE) smaller than 1 cm was accepted for the water level calibration. The bottom friction in the model defined by Nikuradse's equivalent sand roughness (k_s) in the channel is set to 0.1 m. Originally the floodplains are divided into three basic categories: forest, cultivated land and meadows/pastures. In the current study, however, all the floodplain areas are considered to be covered with the same type of vegetation for the sake of simplicity (Fig. 1).

Figure 2 Water level difference along the river channel axis for different floodplain friction values (a), and absolute difference of the flow velocity with and without floodplain friction formulation (b).

The calibration of the floodplain friction with water level measurements can be achieved either with (a) the traditionally used Lindner–Pasche friction law (Eq. 1) in addition to the bottom friction (Fig. 2a, red line) or (b) a higher bottom friction value on the floodplains (Fig. 2a, blue line). The calibration results are shown in Fig. 2. In (a) 5 cm thick stem evenly spaced by 5 m intervals are added to the floodplains with the LIND formulation, while keeping the bottom friction equal to the channel (red line). In (b) k_s=0.5 m is used on the floodplains (blue line). In the case of the same friction value of the river channel (k_s=0.1 m) being used for the floodplains, the momentum is too high and the simulated free surface does not fit the measurements (black line). In both cases the water level RMSE was reduced by a factor of 4 (see values within brackets).

The reader may ask himself or herself about which approach to be used. In this case it is useful to compare the absolute difference of the flow velocity with and without the floodplain friction formulation (see Fig. 2b). It can be seen that while differences in the main channel can be neglected (<0.05 m s⁻¹), those on the floodplains cannot (up to 0.4 m s⁻¹). In other words, when model calibration is based only on measurements in the river channel the hydrodynamics on floodplains is not guaranteed to be correctly simulated. It is important to point out that in the case of unsteady flow conditions a friction formulation dependent on water depth is not desired. However, if flow velocity measurements are not available on the floodplains (usually the case), little can be done in terms of calibration. Furthermore, remote sensing data of vegetation characteristics (light detection and ranging technology) have been used in flood modelling in the last 20 years, but the accuracy of these measurements should also be taken into account (see Cobby et al., 2003; Antonarakis et al., 2008; Dombroski, 2017).

An alternative to the deterministic approach in such situations, when there is a potential shortcoming in the modelling process due to a lack of information, is to carry out an UQ. As explained in Sect. 3, this analysis can be used to estimate the combined effect of all uncertain input parameters that significantly affect model results by means of a probabilistic approach.

4.2 Input parameters

The next step now is to calibrate the remaining floodplain friction formulations with water level measurements. In order to make a comparison to the LIND approach, first H=10.0 m is set to ensure emergent conditions in all formulations. A second scenario is then calibrated for submerged conditions, in which H=1.0 m. Tables 1 and 2 present the calibrated parameters for each one of the scenarios. After the calibration all model results presented RMSE smaller than 1 cm for water level. (The density of stems m is calculated as 1∕d², where d is the distance between stems.)

Table 1 Floodplain friction parameter values calibrated under emergent conditions (H=10.0 m).

Table 2 Floodplain friction parameter values calibrated under submerged conditions (H=1.0 m).

For the UQ it is required that all sensitive parameters relevant to model results should be considered for the determination of the PIs. Once the parameters are chosen a very important step follows: an error or deviation should be carefully assigned to each parameter. This variation should be small enough to be treated as an error, but large enough to include the actual parameter uncertainty (due to the lack of knowledge). Unfortunately there is no general rule for choosing a proper value, since different aspects might contribute, for example, measurement accuracy, spatial/time variances and numerical representation of process. In the current study, the chosen variations for the parameters related to the vegetation species (C_D, LAI, χ) are based on values given in Aberle and Järvelä (2013) and Västilä and Järvelä (2014). For the remaining parameters (H, D, d, K) a standard deviation of 10 % of the calibrated value is assumed (σ=0.1μ). The vegetation height H is only included in the analysis under submerged conditions. Input deviations are treated as errors and, therefore, represented by a Gaussian distribution. This implies that there is a 99.7 % probability that the parameter value is found within $[\mathit{\mu }-\mathrm{3}\mathit{\sigma },\mathit{\mu }+\mathrm{3}\mathit{\sigma }]$.

Table 3 Input parameter deviations and ranges.

^* Only for submerged conditions.

Finally, the UQ methods presented in Sect. 3 are applied with the input uncertainties given in Table 3. The FOSM method is evaluated through central finite difference; hence, 2n model evaluations are necessary (n refers to the number of input variables). For the MC method, a sample size of 1000 was used for the evaluation. Although MC sample sizes are usually considered in the range of 10⁴–10⁵, previous tests with 10⁴ samples showed very little difference in results. The NIPC method (metamodelling) requires fewer samples than MC, because a polynomial function fitted to the samples is then used for the evaluation of results. In this case results with 100 samples for the metamodel were sufficient for approximating MC results.

5 Results and discussion

The numerical model was evaluated with the four floodplain friction formulations. A constant discharge of 7870 m³ s⁻¹ was imposed at the upstream boundary and at the downstream boundary a corresponding free surface based on a discharge curve. Model results were analysed after a steady state was achieved in the simulation and presented in the form of prediction interval (PI) with a 95 % probability of occurrence. It should be noted that the PI represents a range of variation around the mean value, which is not necessarily symmetric (MC and metamodelling). Once Monte Carlo results are available scatterplots can be used to investigate the behaviour of the model. They represent here the projection of the N values of model outputs against the N values of each of the n input parameters. Scatterplots are a very simple and informative form of sensitivity analysis, since they provide a visual representation of the relative importance of the input parameters (Saltelli et al., 2008). In order to summarize the scatterplots with a single value, standardized regression coefficients (SRCs) can be determined through linear regression.

Figure 3 Scatterplots of flow velocity results against floodplain friction input parameters at node 13470.

The scatterplots together with the SRCs are presented in Fig. 3, in which the sampled values of each input parameter (see Table 3) are plotted against the flow velocity results at node 13 470 (see Fig. 1) under emergent conditions. This node is considered representative of the hydrodynamics on vegetated floodplains. The SRC indicates whether there is a significant linear relation between input and output (values given in brackets). For instance the permeability coefficient (K) in BATT (1.0) is linearly related to the flow velocity, due to the approximation used under emergent conditions (see Sect. 2.4). The distance between stems (d) in LIND (0.912) and BAPT (0.805) also presents a significant linear effect on flow velocities. Interestingly the drag coefficient (C_D) presents a much stronger linear relation to the flow velocity in JAER (−0.713) than in BAPT (−0.355). In summary the permeability coefficient (K), the distance between stems (d) and the drag coefficient (C_D) in JAER were identified as inputs from the friction formulations with the largest linear effect on flow velocity. Although the remaining parameters showed no significant linear relation ($\left|\mathrm{SRC}\right|<\mathrm{0.5}$), higher-order dependencies may still exist.

Figure 4 Prediction interval of flow velocity under emergent conditions (rows: friction formulations; columns: UQ methods).

In Fig. 4 the UQ of the flow velocity under emergent vegetation conditions is presented. It can be observed that similar results are obtained with the three UQ methods. Among the friction formulations, LIND and JAER presented smaller variations on average and the PI exceeds 0.2 m s⁻¹ only at the left floodplain in the middle of the river reach. On the other hand the BATT approach appears to be the most sensitive one, followed by BAPT. The approach from Battiato and Rubol (2014) results in PI >0.2 m s⁻¹ on most of the floodplains. Relative to results with calibrated values all formulations showed variations above 10 % on the floodplains. The same analysis is carried out for submerged conditions (see Fig. 5). Because the LIND formulation is only valid for emergent conditions (independent of H), submerged conditions cannot be accounted for. As expected all results present on average larger PI than under emergent conditions, due to the addition of H in the analysis. The floodplain PI of flow velocity is mostly above 0.2 m s⁻¹ and in BATT the PI >0.5 m s⁻¹ is present on floodplains located at the inner bends of the river reach. In the channel the PI in BATT is the largest one and exceeds 0.05 m s⁻¹ all along the upstream river section. Relative to results with calibrated values the variations now exceed 25 % on the floodplains, and in BATT at shallow regions up to 100 %.

Figure 5 Prediction interval of flow velocity under submerged conditions (rows: friction formulations; columns: UQ methods).

As explained in Sect. 3 results from the UQ can also be represented in terms of risk, i.e. a probability of failure (P_f). This is a more suitable analysis for when results must not exceed a given threshold. For instance, a threshold of 0.1 m s⁻¹ above the mean value is used for the analysis of the flow velocity (see Fig. 6). In other words, the probability of exceedance of $\stackrel{\mathrm{‾}}{u}+\mathrm{0.1}$ m s⁻¹ was calculated. Because in the current study the difference among the UQ methods was not significant, results are now presented only from metamodelling. Results indicate that there is a larger probability that velocities are found above $\stackrel{\mathrm{‾}}{u}+\mathrm{0.1}$ with the BATT approach. Under submerged conditions velocities are more likely to exceed this threshold. In BAPT and BATT the probability of failure can be higher than 10 % on the floodplains.

Figure 6 Risk analysis of flow velocity for threshold $\stackrel{\mathrm{‾}}{u}+\mathrm{0.1}$ m s⁻¹ under emergent (a–d) and submerged (e–h) conditions.

Figures 4 and 5 show that results using the UQ methods are similar for this case study. Although the FOSM method is the least expensive alternative among the ones presented (only 2n model evaluations), it should be used with caution as it is a first-order method (linear approximation). In studies involving strongly dominant non-linear processes (e.g. turbulence modelling, sediment transport, unsteady conditions) further comparisons should be carried out. On the other hand, Monte Carlo-based methods have the advantage that the analysis under any conditions is possible. Although a large number of simulations is required for obtaining trustful results, alternatives such as the NIPC make them more feasible by reducing the sample size by at least 1 order of magnitude. For instance, Mouradi et al. (2016) carried out the UQ of a computationally intensive morphodynamic model, to which they applied pure MC and metamodelling methods.

When compared to emergent conditions the overall uncertainty of submerged conditions is significantly larger. This is an expected result in UQ as there is an additional input (vegetation height H) that significantly contributes to model performance. The Lindner–Pasche floodplain friction formulation is by definition only valid for emergent conditions. Thus, a different approach is needed when submerged conditions should be taken into account. Additionally, Brandimarte and Di Baldassare (2012) warn that when simulating flood scenarios attention must be given to parameter compensation between floodplain and channel resistance coefficients, so that reasonable values are chosen.

An important topic not only regarding UQ but numerical simulation in general, is the matter of input uncertainty definition. When performing a numerical simulation that is based on physical processes one will eventually need to validate calculations with measurements. Also, initial and boundary conditions are usually based on measurements of the original process. That is to say one should know a priori how accurate the available measurements are. This is usually not a trivial task, since measurement errors may not be easily evaluated (see Taylor, 1997). For instance, Di Baldassarre and Montanari (2009) published a study that focused only on the uncertainty in river discharge observations. Although it was attributed a standard deviation of 2.7 % for discharge measurement errors, the authors emphasize that this value is associated with their case study, and thus any generalization should be attributed with care. For uncertainties related to floodplain friction there are no such reference studies known to the authors. In that case, a suggested practice is to start with commonly used value ranges in the literature and apply a 6σ range ($[\mathit{\mu }-\mathrm{3}\mathit{\sigma },\mathit{\mu }+\mathrm{3}\mathit{\sigma }]$) for the total parameter variation as a rule of thumb. Of course available experience in the topic of investigation should be also taken into account.

6 Conclusions

A framework for the estimation of uncertainties of hydrodynamic models on floodplains was presented. A traditional resistance formula used for river modelling together with three more recent approaches to floodplain friction were considered for carrying out an uncertainty quantification. The analysis was performed by means of three different methods: traditional MC, FOSM and NIPC (metamodelling). A two-dimensional model of a 10 km reach of the River Rhine was calibrated under steady flow conditions and the analysis was based on flow velocity results.

From the scatterplots it was identified that the permeability coefficient (K), the distance between stems (d) and the drag coefficient (C_D) in JAER produce the largest linear effect on flow velocity among the inputs from the friction formulations. Altogether the tested floodplain friction formulae produced qualitatively similar results, in which uncertainties in flow velocity are most significant where the resistance coefficient was modified. Under emergent conditions, larger velocity variations are obtained with the formulations of BAPT and BATT. Variations from the latter also included the river channel. Under submerged conditions all approaches resulted in larger uncertainties, as the vegetation height was included in the analysis. Although the BATT approach presented once again the largest variations among the analysed methods, results were consistent not only qualitatively, but also quantitatively. In summary, among the tested floodplain friction formulae the JAER approach presented on average the smallest prediction intervals, i.e. the smallest variance. It is important to keep in mind that UQ results are not only dependent on the defined input parameter deviations, but also on the number of parameters considered in the analysis. In that sense, the BATT approach is still attractive as it reduces the current analysis to a single parameter, the canopy permeability K.

The three UQ methods compared gave similar results, which means that FOSM is the most efficient in this case. Despite being a very simple method to apply, FOSM will only produce good results when the first-order approximation is sufficient to describe the sensitivity of the system. In the presented study this was the case, probably because all the chosen inputs are directly correlated to the resistance coefficient. Research on related topics such as floodplain mapping usually focuses on the analysis of uncertainties that relies on Monte Carlo-based methods (e.g. Di Baldassarre et al., 2010; Domeneghetti et al., 2013). Several further topics could be listed here for future development, e.g. unsteady flow, boundary conditions, not to mention sediment transport modelling. However, the most important is first to be aware of the limitations of the available information and tools. Are there enough measurements for an acceptable calibration in the study area? Is the chosen numerical model capable of correctly representing the physical process under the desired conditions? As basic as it may sound, if those questions cannot be answered, any kind of analysis involving uncertainties will fail in providing useful results.