Uncertainty Analysis of a Two-Dimensional Hydraulic Model

: A reliability approach referred to as the point estimate method (PEM) is presented to assess the uncertainty of a two-dimensional hydraulic model. PEM is a special case of numerical quadrature based on orthogonal polynomials, which evaluates the statistical moments of a performance function involving random variables. When applied to hydraulic problems, the variables are the inputs to the hydraulic model, and the ﬁrst and second statistical moments refer to the mean and standard deviation of the model’s output. In providing approximate estimates of the uncertainty, PEM appears considerably simpler and requires less information and fewer runs than standard Monte Carlo methods. An example of uncertainty assessment is shown for simulated water depths in a 46 km reach of the Richelieu River, Canada. The 2D hydraulic model, H2D2, was used to solve the shallow water equations. Standard deviations around the mean water depths were estimated by considering the uncertainties of three main input variables: ﬂow rate, Manning’s coefﬁcient and topography. Results indicate that the mean standard deviation is <27 cm for the considered ﬂow rates of 759, 824, 936, 1113 m 3 /s. Higher standard deviations were obtained upstream of the topographic shoal at the municipality of Saint-Jean-sur-Richelieu. The PEM method adds further value to the H2D2 model predictions as it indicates the magnitude and the spatial variation in uncertainties. The effort required to complete an uncertainty analysis using the PEM method is minimal and the resulting insight is meaningful. This knowledge should be incorporated into decision-making in the context of ﬂood risk assessment. wind forces and the Coriolis effect were not considered. This is because there was insufficient information to describe the uncertainty in the wind data available for the Richelieu River. The resulting water depths were calculated by subtracting water surface elevations from river bed elevations. The model was used to evaluate the water levels and flow in the Richelieu River. The water level information could then be utilized by other studies relating to the impact of water levels in the river. In all of these applications no quantification of uncertainty in the model outputs was made. The model took approximately one hour to provide a steady-state solution on a computer with a core i7 processor of 3.6 GHz and 16 Gigabytes of RAM (Random access memory).


Introduction
Flood inundation models numerically transform discharge rates into predictions of water depth, flow velocity, etc. [1]. The products of hydrodynamic models are often taken as inputs into other models studying pollutant or sediment transport. Also, studies have applied such models in both rural and urban environments [2,3]. Although deterministic, they involve numerous uncertainties, mainly from the input variables, which are due to the cumulative nature of the models and can lead to considerable uncertainty in the results. The origin of the uncertainty is in the simplifying nature of models and our inability to precisely define the required input quantities for the models. Knowledge of the type and magnitude of the uncertainties is crucial to understand and interpret a model's results. Therefore, uncertainty analysis, conducted to quantify the uncertainty in the model's outputs values induced by the uncertainty in the inputs, increases the confidence in hydraulic model predictions [4].
It indicates the magnitude and spatial variation of the uncertainty associated with model results, and helps hydrologists and decision makers to identify the needed improvements. Uncertainty in flood inundation models has been addressed in several studies, focusing on reviewing uncertainties in the input data, model parameters and model structure [5][6][7][8][9][10]. The uncertainties in the model structure are often neglected as there is no general agreement on the best approach to consider model uncertainty. On the other hand, uncertainties in the input data and model parameters are generally taken into consideration, as is the case in this study. There are a number of simplified approaches to flow routing that have been applied to flood modelling, as well as to the contaminant transport problem [11][12][13][14][15][16].
Uncertainty analysis provides additional information to the computations performed by the model. The value of the hydrodynamic modelling outputs is increased by the description of the uncertainty in those outputs. It is possible to place an uncertainty measure (standard deviation or confidence interval) around the computed water depth. This information adds value to the hydrodynamic modelling output and summarizes all of the knowledge that the engineer has on the model output. Uncertainty analysis may help to improve model inputs and model parameters before it is applied the problem. It can be used to complement model calibration by providing insights into how good the model actually is. Better decisions can be made in hydrology and water resource management by taking into account the uncertainties in hydrodynamic modeling.
The source of the uncertainty can be attributed to different aspects of the modelling process. In this study, the contribution of the three most commonly considered potential sources of uncertainty by flood inundation models is analyzed: discharge rate, Manning's coefficient and topography. Each of these input variables and parameters is subject to measurement, modelling errors and/or sampling errors [17,18]. Uncertainty in the flow rate originates from the stage-discharge relation and the flood frequency analysis. The estimation of flood frequencies is well known as a key task in flood hazard assessment. Flood frequency analysis is based on the presumption that the observed flood discharges come from a parent distribution and can be described by a probability distribution. Water levels are a crucial type of data required for hydrodynamic modelling. Water levels are used as boundary conditions and also as observation data in the calibration of a model. The uncertainty associated with the measured water level is not evaluated in the uncertainty analysis because it is measured with very high accuracy. Based on the literature, the errors reported in the literature for water levels have an accuracy of 1 cm [19]. To generate the probability of a given flood level and its respective flow rate, the flood frequency analysis focuses on the statistical analysis of peak discharge by fitting a statistical model directly to the observed sample of peak flow data [1]. The uncertainty in a T-year event discharge is estimated using the confidence intervals. The T-year event is assumed to be normally distributed, which is asymptotically true for most quantile estimators [20]. The uncertainty is then expressed with the following confidence interval: x T − Z 1−p/2 var(x T ); x T + Z 1−p/2 var(x T ), (1) where Z p is p-th quantile of the standard normal distribution and var(x T ) is the variance of the T-year event x T . Some studies have attempted to quantify boundary condition uncertainty in hydraulic modelling [21]. Manning's coefficient is usually used as the calibration parameter in hydraulic modelling [1,22]. This coefficient represents the resistance to flood flows in channels and floodplains. Various studies point out that hydraulic models can be very sensitive to these coefficients [23]. As it cannot be measured, it is estimated empirically or indirectly by laboratory and field methods, which necessarily imply a significant degree of uncertainty [24,25]. Remotely sensed data may provide ways of estimating physical components of the friction term [24]. The variance in Manning's n was estimated using a reference that investigated the uncertainty in hydraulic parameters [26]. Based on experiments and field observations, different error probability distributions and variation coefficients are suggested in the literature to characterize Manning's coefficient, such as normal [27][28][29], lognormal [30,31], or uniform distributions [26,32]. The authors concluded that coefficients of variation for Manning's n can range from 0.1 to 0.3.
Knowledge of the surface topography is probably the most important for the accuracy of hydraulic models [33]. Lack of adequate topographic data can lead to problems for the description of flooded areas provided by the hydraulic model [8,34,35]. Riverbanks and floodplains are commonly identified by digital elevation models (DEM), cross section, or geometrical description of hydraulic structures [19,36]. DEM provides continuous elevation information in 2D models, including floodplain morphology and river bathymetry [37]. The errors in the DEM can span a few orders of magnitude, e.g., low precision DEMs show standard deviations ranging from 5 m to 50 m, whereas for high precision LiDAR data (light detection and ranging) the error is decreased between 0.15 m and 0.3 m [38]. The elevation errors are generally assumed to follow a normal distribution [39][40][41]. The topographic uncertainty in the study [33] used only vertical accuracy as the source of uncertainty.
Other topographically-derived model attributes such as geometry and spatial resolution can also have significant impact on the overall hydraulic model output [42][43][44][45]. Hunter and Bates [24] pointed out that modern LiDAR DEMs are sufficiently accurate and resolved for simulating floods in urban areas.
The uncertainty of any model output (water surface elevations and inundation extent) is obtained by combining the standard deviation of the following variables: discharge rate, Manning's coefficient and topography [46,47]. Numerous approaches exist to quantify the resulting model uncertainties [48,49], including analytical methods [50], generalized likelihood uncertainty estimation (GLUE) [51,52], parameter estimation [53], approximation methods (e.g., first-order second-moment method), Melching [54], simulation and sampling-based methods [55], Bayesian methods [56], analytical methods, methods based on the analysis of model errors, methods based on fuzzy sets [57,58], etc. The selection of a suitable method depends on the knowledge of the river's hydraulics, the availability of data, and the modelling approach, as suggested in the code of practice by Pappenberger and Beven [59]. Hydrodynamic models involve differential equations that are too complex to solve analytically and therefore analytical uncertainty analysis methods are not applicable. When using complex models or functions, approximate uncertainty analysis methods are practical.
Until now, in open channel hydraulics, uncertainty analysis has been limited to models that are not as complex as two-dimensional hydrodynamic models [60]. The majority of the approaches for uncertainty analysis are stochastic and most commonly based on the Monte Carlo method. However, this type of analysis requires significant computer processing time, which can become impractical for complex models like a two dimensional hydrodynamic model [22,61]. The simulations solve the more complex equations of mass continuity and conservation of momentum in two-dimensions using numerical methods. The data requirements are more extensive and have greater potential to create uncertainty in the output. There are methods capable of quantifying the uncertainty in two-dimensional hydrodynamic models with cheap, efficient computing power readily available. In such cases, the point estimate method (PEM), which gives approximate results of the uncertainty with significantly less computational effort, seems to be a promising alternative. PEM represents a special case of numerical quadrature based on orthogonal polynomials, allowing the evaluation of the low order of performance functions of independent random variables. For normal variables, the point estimate method corresponds to Gauss-Hermite quadrature. The difficulty is that PEM does not provide the tails of the output variable's distribution, whereas the Monte Carlo analyses usually generate continuous probability distribution functions [26,62]. The probability distribution functions for each of the model inputs were specified by assuming that the uncertainty in each model input is normally distributed. Normal distribution has been utilized in other similar studies [7].
There has not been research undertaken to quantify the uncertainty in hydrodynamic modelling outputs and to identify the sources of the uncertainty. It is favorable to the decision maker to have the uncertainty in the model output quantified. The objective of this paper is to demonstrate the applicability of the PEM for computing of the mean and standard deviation of the simulated water depths. Intrinsic uncertainties of three main input variables were considered: flow rate, Manning's coefficient and topography, for four different flow rates. An uncertainty assessment was conducted for a 46 km long reach of the Richelieu River, Canada. The finite-element 2D hydraulic model, H2D2, was used to solve the shallow water equations. A discussion of how the uncertainty results can be used is provided at the end.

Methods
The uncertainty analysis method proposed involved the following steps: (i) define the probability distribution for each of the input variables: X 1 (flow rate), X 2 (Manning's coefficient) and X 3 (topography); (ii) generate sampled inputs for the identified variables by using PEM; (iii) run the hydraulic model based on these input parameters (recalculate the model for each sampled input); and (iv) quantify the uncertainty of the model outputs by computing the mean and standard deviation of the water depths (Table 1). was used to solve the shallow water equations. A discussion of how the uncertainty results can be used is provided at the end.

Methods
The uncertainty analysis method proposed involved the following steps: (i) define the probability distribution for each of the input variables: X (flow rate), X (Manning's coefficient) and X (topography); (ii) generate sampled inputs for the identified variables by using PEM; (iii) run the hydraulic model based on these input parameters (recalculate the model for each sampled input); and (iv) quantify the uncertainty of the model outputs by computing the mean and standard deviation of the water depths (Table 1).

Preprocessing Processing Postprocessing
PEM approach sampling H2D2 model Propagation of uncertainty

Study Area
The source of the Richelieu River is Lake Champlain in the USA and it empties northward into the St. Lawrence River, Canada. The discharge at which Lake Champlain drains at increased rates is controlled mainly by an area of rocky shoals located at St-Jean-sur-Richelieu, about 36 km downstream ( Figure 1). A constriction in the Richelieu regulates the outflow from Lake Champlain and dictates the upstream water level variations. The flow contributions from tributaries were ignored during the calibration process because of the significant storage capacity of Lake Champlain, which contributes to more than 95% of the Richelieu River's flow into the St. Lawrence River [63].

Study Area
The source of the Richelieu River is Lake Champlain in the USA and it empties northward into the St. Lawrence River, Canada. The discharge at which Lake Champlain drains at increased rates is controlled mainly by an area of rocky shoals located at St-Jean-sur-Richelieu, about 36 km downstream ( Figure 1). A constriction in the Richelieu regulates the outflow from Lake Champlain and dictates the upstream water level variations. The flow contributions from tributaries were ignored during the calibration process because of the significant storage capacity of Lake Champlain, which contributes to more than 95% of the Richelieu River's flow into the St. Lawrence River [63].

Hydraulic Model of the Richelieu River
The simulations are based on a two-dimensional finite element hydraulic model of the Richelieu River using the H2D2 software [64]. H2D2 solves the 2D depth-averaged shallow-water equations (SWE) based on the conservative form of mass Equation (2) and the two equations of momentum conservation (3) and (4) under steady-state flow conditions with special treatment of drying-wetting areas.

Hydraulic Model of the Richelieu River
The simulations are based on a two-dimensional finite element hydraulic model of the Richelieu River using the H2D2 software [64]. H2D2 solves the 2D depth-averaged shallow-water equations (SWE) based on the conservative form of mass Equation (2) and the two equations of momentum conservation (3) and (4) under steady-state flow conditions with special treatment of drying-wetting areas.

Hydraulic Model of the Richelieu River
The simulations are based on a two-dimensional finite element hydraulic model of the Richelieu River using the H2D2 software [64]. H2D2 solves the 2D depth-averaged shallow-water equations (SWE) based on the conservative form of mass Equation (2) and the two equations of momentum conservation (3) and (4) under steady-state flow conditions with special treatment of drying-wetting areas.
∂ ∂x ∂ ∂x where q x and q x are the flow rates in x and y direction (m 3 /s), h is the water level (m), H is the depth of water, c is wave velocity (m/s), r is the water density, ρ is the density, u(u, v) are the components of the of velocity vector (m/s), f c is the Coriolis force, τ ij is the Reynolds stresses (m/s 2 m), τ b x and τ b y are the bottom frictions in the x and y directions (kg/s 2 m), τ s x and τ s y are the surface frictions in the x and y directions (kg/s 2 m), and x(x, y) are the components of coordinate orientation of X(m). The model does not take into account the influence of the wind and the Coriolis force. The continuous fields h, u and v are discretized over a structured mesh of finite element triangles of six nodes, referred to as T6L.
The H2D2 model is based on the assumptions of incompressibility, hydrostatic pressure, and stable riverbed. The finite-element discretization has been detailed by [65,66]. The H2D2 model calculates the flow velocity, discharge and water levels in a 2D domain using an implicit time scheme and direct Newton Raphson steps. In doing so, it also takes into account the wetting and drying of the river banks, based on the water flow and level [67]. The finite element mesh stores all of the input variables needed for the resolution of St-Venant equations, as well as the resulting variables for the 2D flow simulation (water level and velocities). The advantage of using this software is its simplicity in terms of data requirements and most particularly the possibility of modifying the river geometry, qualities that have been demonstrated in several flood inundation studies [68][69][70].
First, H2D2 was applied in the Richelieu River watershed located about 20 km east of Montreal, Canada ( Figure 1). The model had been pre-calibrated by Environment Canada with the high water event of 6 May 2011, under steady-state flow conditions without wind forcing [71]. The flow during the survey was at an approximate rate of 1550 m 3 /s at the Fryers Rapids gauging station. Calibration was carried out by adjusting the Manning's coefficient. There was a difference between simulated and measured water levels at Saint-Jean-sur-Richelieu of 1cm. The results for the calibration for this high event on 6 May 2011 are presented in Figure 3.  , and x(x, y) are the components of coordinate orientation of X(m). The model does not take into account the influence of the wind and the Coriolis force. The continuous fields h, u and v are discretized over a structured mesh of finite element triangles of six nodes, referred to as T6L.
The H2D2 model is based on the assumptions of incompressibility, hydrostatic pressure, and stable riverbed. The finite-element discretization has been detailed by [65,66]. The H2D2 model calculates the flow velocity, discharge and water levels in a 2D domain using an implicit time scheme and direct Newton Raphson steps. In doing so, it also takes into account the wetting and drying of the river banks, based on the water flow and level [67]. The finite element mesh stores all of the input variables needed for the resolution of St-Venant equations, as well as the resulting variables for the 2D flow simulation (water level and velocities). The advantage of using this software is its simplicity in terms of data requirements and most particularly the possibility of modifying the river geometry, qualities that have been demonstrated in several flood inundation studies [68][69][70].
First, H2D2 was applied in the Richelieu River watershed located about 20 km east of Montreal, Canada ( Figure 1). The model had been pre-calibrated by Environment Canada with the high water event of 6 May 2011, under steady-state flow conditions without wind forcing [71]. The flow during the survey was at an approximate rate of 1550 m s ⁄ at the Fryers Rapids gauging station. Calibration was carried out by adjusting the Manning's coefficient. There was a difference between simulated and measured water levels at Saint-Jean-sur-Richelieu of 1cm. The results for the calibration for this high event on 6 May 2011 are presented in Figure 3. The influence of the wind forces and the Coriolis effect were not considered. This is because there was insufficient information to describe the uncertainty in the wind data available for the Richelieu River. The resulting water depths were calculated by subtracting water surface elevations from river bed elevations. The model was used to evaluate the water levels and flow in the Richelieu River. The water level information could then be utilized by other studies relating to the impact of water levels in the river. In all of these applications no quantification of uncertainty in the model outputs was made. The model took approximately one hour to provide a steady-state solution on a computer with a core i7 processor of 3.6 GHz and 16 Gigabytes of RAM (Random access memory). The influence of the wind forces and the Coriolis effect were not considered. This is because there was insufficient information to describe the uncertainty in the wind data available for the Richelieu River. The resulting water depths were calculated by subtracting water surface elevations from river bed elevations. The model was used to evaluate the water levels and flow in the Richelieu River. The water level information could then be utilized by other studies relating to the impact of water levels in the river. In all of these applications no quantification of uncertainty in the model outputs was made. The model took approximately one hour to provide a steady-state solution on a computer with a core i7 processor of 3.6 GHz and 16 Gigabytes of RAM (Random access memory).
The topography data for the model was determined using a DEM with a cell size of 1 × 1 m 2 and an accuracy of ±0.3 m derived from LiDAR datasets obtained in 2008, 2010, and 2011, provided by Environment Canada (EC). Based on the various classes of point from the data sources, EC used only the ground class, which represents 85% of the LiDAR points surveyed. The other classes, such as the building and vegetation, were not used in the DEM. A set of orthophotos from the 2009 database for the administrative region of Montérégie, Quebec was used to validate the limits of the water courses. The elevations of the points below the water surface came from bathymetry data. Bathymetric information describes the geometry of a river. The bathymetry of the riverbed was taken from the Canadian Hydrographic Service (low river) and bathymetric transects were provided by Parks Canada (mid river). To import topobathymetry into the model, the individual soundings were interpolated onto the finite element grid nodes. The interpolation of each individual sounding reduces the amount of uncertainty in the bathymetry in the model. Uncertainties in the bathymetric data were not explicitly dealt with in this study. The range of calibrated Manning's coefficients between of 0.02-0.036 was used as the mean value of the random field over the study area. The substratum upstream consists of fine sand and silt. The substrate is much coarser with some boulders, approaching the shoal in Saint-Jean-sur-Richelieu.
The upstream boundary condition of the model is the flow rate at Rouses Point and stage conditions used at the downstream limit of the model near the Fryers Rapids station (EC 02OJ007). The water level was estimated using the established stage-discharge relationships. Q and H are, respectively, the discharge (m 3 /s) and water level (m)

Point Estimate Method
Knowledge of the mean and the standard deviation of each input variable at each of the model's grid cells is necessary to compute the combined uncertainty of the outputs using the PEM. The method estimates the mean and standard deviation of a performance function using statistical moments of the random input variables. Originally proposed by Rosenblueth [72], the method has been improved in the past decade [73,74]. The objective is to evaluate the model at a discrete set of points in the uncertain parameter space and to combine uncertainties of the random variables with proper weighting in order to calculate the mean and standard deviation of the output, which is the water depth in this case. These outputs provide useful information for researchers and managers, but in the absence of uncertainty information, the mean values alone do not show the expected range of possible water depth values.
The standard PEM considers a set of random variables X i (X 1 , X 2 . . . , X n ) with a probability distribution function f Xi (x), and a performance function Y, obtained from the deterministic relationship of X (Y = g(X i )) to represent the output variable Yi. In this application, the considered independent random variables include the flow rate as the upstream boundary condition, Manning's coefficient, and topography. In this case, the performance function is the water depth, g(x). The input variables are stochastic, i.e., they can be described in terms of their mean and variance. The PEM method is used to approximate the mean, variance, and/or higher moments of the outputs. The accuracy of the prediction of the uncertainty in the model outputs depends largely on the ability to estimate the uncertainty in the inputs.
The statistical moments of the model output can be calculated in terms of the respective moments of the input as follows: where u xm is moment of order m and µ X is the first moment of f X (x). The integrals given in Equation (4) are solved by applying the approximate integration with PEM. The continuous random variable X is replaced by a discrete random variable with a probability mass function (pmf) of the same moment order m as f x (x). The pmf p X (x) is transformed using g(X) to obtain another discrete function with pmf p Y (y). In the continuous case, p Y (y) is used to calculate the moments of Y. In the case of a 1D integral, the mth central moment of f X (x), µ Xm , is obtained as follows: The optimal coordinates are those at which the integral and the corresponding weights are selected by the Gaussian quadrature procedures. In fact, the PEM can be considered as an application of the Gaussian quadrature procedures [75]. Table 2 shows the coordinates along the abscissa and weights for the three first quadrature formulas based on the normal distribution as the weight function. For a case in which the performance function involves n independent random variables assumed to be approximately Gaussian, Rosenblueth suggested that x can be estimated at more than two points. The three-point method applied herein considers a central point at x = µ x and two points x + and x − symmetrically distributed around the mean. Designating the weight of the central point as P and of the other two points as P_ and P+, the respective weights can be written as: The solutions of these equations gives the following results: In most cases, the calculations are made at two points and the expected value of the k-th power of the probability distribution of Y is determined by: where Y is a deterministic function of X = g(X); E[Y m ] is the expected value of Y to the power m; y + is the value of Y at the upper point x + , which is greater than the mean µ x ; y − is the value of Y at the lower point x − , which is less than the mean µ x ; y m is the value of Y at the mean µ x . If m is large, the method requires a considerable computational effort, as is the case with 2D SWE models.
To assess the uncertainty propagation, the performance of the model is evaluated through (n + 1) p simulations, where: p is the number of uncertain variables and n is number of nodes. Each simulation corresponds to the perturbation of a given random variable with a ± √ 3 standard deviation about its mean value. The point estimate method needed only 27 calculations of the H2D2 model for each flow case. The PEM method is relatively simple in formulation when compared to other uncertainty analysis methods, such as the Monte Carlo Method. The disadvantage of Monte Carlo analysis is that it can be computationally demanding. This can be an issue for complex models that require a significant amount of time to compute. Figure 4 shows the procedure applied to estimate the statistical moments of the model output for three input variables. For the purpose of numerical computation, the random fields of the uncertain parameters are represented in terms of uncorrelated random variables. The H2D2 model of the Richelieu River calculates values for Equation (13) at all the nodes within the model domain. PEM simulations were conducted by using three variables with the proper weighting ( Table 3). The sum of the all weights is equal to 1. The disadvantage of PEM is that it loses precision as the nonlinearity of the model increases and if the moments higher than the second moment are to be obtained, and the fast growing of the number of calculations needed if the number of random variables increases, which increases the complexity of the analyses and the computation time [74]. This is especially the case if the input parameters are correlated, requiring the computation of covariance for each pair of variables and inclusion of the co-variances in the uncertainty calculations Therefore, PEM gave reliable and accurate results with the third lowest moments. It is relatively easy to run with the bulk of the effort required to utilize the method. The disadvantage of the PEM method, is that it does not provide a measure of the influence of each random variable to the overall variance, so it is not an adequate method to filter the most relevant random variables. Also, the PEM method does not take into account the distribution of the input parameters. It relies only on the mean and the variance to describe the parameter. Some of the shortcomings of the method arise from the reality that the PEM method is an approximate uncertainty analysis method.  Table 3). The sum of the all weights is equal to 1. The disadvantage of PEM is that it loses precision as the nonlinearity of the model increases and if the moments higher than the second moment are to be obtained, and the fast growing of the number of calculations needed if the number of random variables increases, which increases the complexity of the analyses and the computation time [74]. This is especially the case if the input parameters are correlated, requiring the computation of covariance for each pair of variables and inclusion of the covariances in the uncertainty calculations Therefore, PEM gave reliable and accurate results with the third lowest moments. It is relatively easy to run with the bulk of the effort required to utilize the method. The disadvantage of the PEM method, is that it does not provide a measure of the influence of each random variable to the overall variance, so it is not an adequate method to filter the most relevant random variables. Also, the PEM method does not take into account the distribution of the input parameters. It relies only on the mean and the variance to describe the parameter. Some of the shortcomings of the method arise from the reality that the PEM method is an approximate uncertainty analysis method.

Flow Rate
The reference flow rates for the hydraulic model were determined by analyzing the frequency of daily recordings at the Fryers Rapids station (between 1972 to 2011). To fully investigate the uncertainty in potential flood events, four flow rates were imposed to constrain the model: 759, 824, 936, 1113 m 3 /s. These values are related to four return periods: 1.25, 1.4, 2, and 5 years used by Quebec's public safety ministry. The hydrologic uncertainty was determined with the standard error of the design flow rate, assuming that its distribution is normal [20] (Table 4). The perturbed discharge Q t is calculated as the mean discharge Q base ± √ 3 standard deviation.
The set of three values for each flow rate Equation (14) was used to obtain the perturbed flow rate considered as an upstream boundary condition for assessment of the performance function (Table 4).

Manning's n Coefficient
The basic roughness parameters were derived from tabulated values [76] and further modified during the calibration the model. The Manning's n coefficient is assumed to be normally distributed about the mean values calibrated during the calibration process of the model [26][27][28] and not through field measurements or observations. The mean n values vary spatially across the domain within the range of 0.02 to 0.036, inferred from substrate samples, and the standard deviation was calculated based on the 3σ rule as 0.0026 (s/m 1/3 ) [77].
The uncertainty in the Manning's n is likely smaller. The reason for this is the calibration process, which makes use of river water levels to indirectly determine the appropriate Manning's n value. Still, the calibration process is not exact and even though these parameters are optimally determined through calibration, they could still be incorrect or have some degree of uncertainty in them. There are alternative parametrization methods that perform almost as well as the optimal Manning's n values. There is also the possibility that the Manning's n parameters may be compensating for errors in other aspects of the model, such as channel bathymetry.
The three random fields generated by the input of the Manning's coefficient in the H2D2, used to evaluate the model's performance, are provided Figure 5. It involves changing the Manning's coefficient along the river by the √ 3 error. For example, this change happens in the area with sills or where the bathymetry suggested a coarser substrate. An increase in the coefficient (upper scenario) yields a lower flow and a higher water level. However, the decrease in the coefficient (lower scenario) includes shallow water depth, faster velocities, and supercritical flows.

Topography
To investigate the effect of topography uncertainty on model predictions, a standard deviation related to each grid cell was generated by kriging the DEM used as the mean value. As was the case in the evaluation of the Manning's coefficient, it was hypothesized that the errors related to the estimation of topography follow a normal distribution and are independent [39][40][41]. Kriging was performed by applying a locally adaptive exponential variogram to the entire study area with a sill value of 0.8 m and nugget effect of 0.03 m (0.18 = 0.03). The sill implies the variance of the data. (b) nominal scenario (µ); and (c) higher scenario µ + √ 3σ .

Topography
To investigate the effect of topography uncertainty on model predictions, a standard deviation related to each grid cell was generated by kriging the DEM used as the mean value. As was the case in the evaluation of the Manning's coefficient, it was hypothesized that the errors related to the estimation of topography follow a normal distribution and are independent [39][40][41]. Kriging was performed by applying a locally adaptive exponential variogram to the entire study area with a sill value of 0.8 m 2 and nugget effect of 0.03 m 2 0.18 2 = 0.03 . The sill implies the variance of the data. The nugget could represent the measurement error or some variation at small scale. The error represents the vertical accuracy DEM (±15 cm), and the vertical error specific to the airborne positioning system (3 cm). The kriging variances were then transformed into standard deviation (i.e., the square root of estimation variance). Figure 6a shows the highest standard deviation with a maximum value of 0.21 m. These changes are incorporated in the steady-state hydraulic model. To quantify the topography's uncertainty, the kriged standard deviation was added to and subtracted from the DEM (mean topography) to obtain the upper and lower confidence levels (± √ 3 standard deviation). With this change, the topography behaves as a flexible body because it will move up or down from a normal distribution and therefore lead to the change in the simulated water surface elevations from the hydraulic modelling. The random topography fields generated by perturbing the mean topography by the √ 3 error are given in Figure 6b-d. With these perturbations, the topography was moved up and down from the mean values, which resulted in a change of the simulated water surfaces (without varying the flow rate and Manning's coefficient in the hydraulic model).

Uncertainty in the Model Output
Having defined the individual uncertainty in the three input variables, PEM simulations were performed in combination with H2D2 to assess the uncertainty which is propagated into the output values. For a given flow rate, the hydraulic model was only run 27 times, in agreement with the 3 × 3 × 3 combinations associated with the μ and μ ± √3σ random fields of the three input variables. The mean and the standard deviation of the water depth for each grid cell were calculated with Equation

Uncertainty in the Model Output
Having defined the individual uncertainty in the three input variables, PEM simulations were performed in combination with H2D2 to assess the uncertainty which is propagated into the output values. For a given flow rate, the hydraulic model was only run 27 times, in agreement with the 3 × 3 × 3 combinations associated with the µ and µ ± √ 3σ random fields of the three input variables. The mean and the standard deviation of the water depth for each grid cell were calculated with Equation (11). Therefore, to fully investigate the uncertainty distributions, 108 model runs (4 × 27) were sufficient to obtain the expected mean and standard deviation of the water depth for the four reference flow rates 759, 824, 936, 1113 m 3 /s. The H2D2 model was prepared to explore the transfer of variability from the Manning's coefficient, flow rate and topography to water depth under steady-state application, without non-linear perturbations. Each run of the H2D2 model had an average time duration of 1 h, which makes flood uncertainty analysis with Monte Carlo unfeasible from a practical point of view in engineering. Still, an approximate uncertainty analysis can be performed with the help of the point estimate method. The mathematical operations defined in Equation (13) were performed within the GIS using the simulated layers with H2D2. The resulting spatial distributions of the mean and the standard deviation of the outputs are shown in Figures 7 and 8, respectively, which were obtained by implementing the point estimate calculations within the GIS software. The expected value or mean of the water depths was determined by running the model with the mean input parameter values.    It can be observed in Figures 7 and 8 that the expected water depths were below 12.06 m deep, whereas the overall standard deviations were less than 27 cm for the four flow rates. The highest water depths were located downstream of Saint-Jean-sur-Richelieu. At Saint-Jean-sur-Richelieu, the river becomes narrower and the bottom rises to the head of the rapids near the shoal. In this reach the slope of the river is quite steep where the riverbed is composed of a coarser substrate. Below the shoal, the slope of the river becomes flatter until the downstream limit is reached. The upstream reach, from Rouses Point to Saint-Jean-sur-Richelieu, is characterized by lower water depths. The river is wide and the impediment to flow is not great, it is a lake-like section of the river with a very minor gradient. This region can be considered as an extension of Lake Champlain. The spatial distribution of the standard deviations shows higher standard deviations in the upstream close to the shoal in Saint-Jean-sur-Richelieu, and lower standard deviations in the far upstream (Figures 8). The increase in uncertainty in water depth is directly proportional to the increase in uncertainty in these model parameters. This finding is expected due to the equation used for calculating the uncertainty, Equation (12). Water depth upstream of Saint Jean is therefore influenced almost entirely by the shoal. The lower standard deviation near the model boundaries is probably in part due to the fact that the water level is controlled at these locations directly at the downstream boundary and indirectly by the It can be observed in Figures 7 and 8 that the expected water depths were below 12.06 m deep, whereas the overall standard deviations were less than 27 cm for the four flow rates. The highest water depths were located downstream of Saint-Jean-sur-Richelieu. At Saint-Jean-sur-Richelieu, the river becomes narrower and the bottom rises to the head of the rapids near the shoal. In this reach the slope of the river is quite steep where the riverbed is composed of a coarser substrate. Below the shoal, the slope of the river becomes flatter until the downstream limit is reached. The upstream reach, from Rouses Point to Saint-Jean-sur-Richelieu, is characterized by lower water depths. The river is wide and the impediment to flow is not great, it is a lake-like section of the river with a very minor gradient. This region can be considered as an extension of Lake Champlain. The spatial distribution of the standard deviations shows higher standard deviations in the upstream close to the shoal in Saint-Jean-sur-Richelieu, and lower standard deviations in the far upstream ( Figure 8). The increase in uncertainty in water depth is directly proportional to the increase in uncertainty in these model parameters. This finding is expected due to the equation used for calculating the uncertainty, Equation (12). Water depth upstream of Saint Jean is therefore influenced almost entirely by the shoal. The lower standard deviation near the model boundaries is probably in part due to the fact that the water level is controlled at these locations directly at the downstream boundary and indirectly by the influence of the upstream boundary defined by the constant flow rate. It can also be observed that an increase of the flow rate from 759 m 3 /s (Figure 8a) to 1113 m 3 /s (Figure 8d) slightly decreases the standard deviations of the water depth values. This can be explained by the hydraulic model for the Richelieu River, which was solely calibrated for the high-water event of 6 May 2011, using a flow rate of 1550 m 3 /s.
These findings indicate that the water depths upstream of Saint-Jean-sur-Richelieu are therefore highly sensitive to the shoal's topography at Saint-Jean-sur-Richelieu. The upstream water slope therefore remains a function of the streambed roughness, since the rivers primary level and rate of discharge is limited to its capacity at the shoal. From the point of view of the engineer that has to build a model, the location of these zones with higher variability can be useful to make decisions such as where to direct efforts to efficiently improve the model. The higher uncertainty located at the reach upstream may be explained by the presence of the Saint-Jean-sur-Richelieu shoal. The latter contains some anthropogenic submerged structures such as steel traps for fishing nets and mill races (Figure 9). At low flow rates, these structures fill with water without contributing to the flow because the velocities are very low, and probably act as a natural sills at low flow rates. Therefore, they are a significant contributing factor to the hydraulic response of this basin. The submerged structures would require relocation because of increased water depths and generally unsuitable conditions for trapping created by dredging and structure operations. The important effect of the shoal on low flow regimes was the reason that the hydraulic model for the Richelieu River was calibrated for the high water event of 6 May 2011. For this reason, and in order to build an accurate 2D hydraulic model, the topography of the model should be carefully considered, especially in the shoal, possibly by acquiring a new river bathymetry. influence of the upstream boundary defined by the constant flow rate. It can also be observed that an increase of the flow rate from 759 m 3 /s (Figure 8a) to 1113 m 3 /s (Figure 8d) slightly decreases the standard deviations of the water depth values. This can be explained by the hydraulic model for the Richelieu River, which was solely calibrated for the high-water event of 6 May 2011, using a flow rate of 1550 m 3 /s. These findings indicate that the water depths upstream of Saint-Jean-sur-Richelieu are therefore highly sensitive to the shoal's topography at Saint-Jean-sur-Richelieu. The upstream water slope therefore remains a function of the streambed roughness, since the rivers primary level and rate of discharge is limited to its capacity at the shoal. From the point of view of the engineer that has to build a model, the location of these zones with higher variability can be useful to make decisions such as where to direct efforts to efficiently improve the model. The higher uncertainty located at the reach upstream may be explained by the presence of the Saint-Jean-sur-Richelieu shoal. The latter contains some anthropogenic submerged structures such as steel traps for fishing nets and mill races ( Figure  9). At low flow rates, these structures fill with water without contributing to the flow because the velocities are very low, and probably act as a natural sills at low flow rates. Therefore, they are a significant contributing factor to the hydraulic response of this basin. The submerged structures would require relocation because of increased water depths and generally unsuitable conditions for trapping created by dredging and structure operations. The important effect of the shoal on low flow regimes was the reason that the hydraulic model for the Richelieu River was calibrated for the high water event of 6 May 2011. For this reason, and in order to build an accurate 2D hydraulic model, the topography of the model should be carefully considered, especially in the shoal, possibly by acquiring a new river bathymetry.

Conclusions
Water depth uncertainties were quantified as the output of the two-dimensional hydraulic model. The uncertainty propagated from the input variables to the model's output was calculated using the point estimate method (PEM). To demonstrate the practical applicability of the PEM for uncertainty analyses, the hydraulic model H2D2 was applied over a 46 km long reach of the Richelieu River, Canada. Three random variables were considered to be sources of uncertainty: flow rate, Manning's coefficient and topography. The combined uncertainty was calculated using the PEM analysis method.

Conclusions
Water depth uncertainties were quantified as the output of the two-dimensional hydraulic model. The uncertainty propagated from the input variables to the model's output was calculated using the point estimate method (PEM). To demonstrate the practical applicability of the PEM for uncertainty analyses, the hydraulic model H2D2 was applied over a 46 km long reach of the Richelieu River, Canada. Three random variables were considered to be sources of uncertainty: flow rate, Manning's coefficient and topography. The combined uncertainty was calculated using the PEM analysis method.
Reasonable estimations of the mean and standard deviation values of water depths were obtained with much less effort using PEM method. Also, the results revealed a standard deviation ranging from 0 cm to 27 cm for the Richelieu River. Regarding the maximum standard deviation, a maximum uncertainty of 27 cm was found mostly upstream of the Richelieu River near the shoal. This is an indication that the model needs to be improved in this area. The results presented showed that the PEM method can be applied with the H2D2 model to perform uncertainty analysis. The PEM method requires less information to describe the model input, fewer model executions, and the resulting insight is meaningful.
Independently of the method used to determine the output's uncertainty, the uncertainty of the model's inputs must be accurately estimated. Although this study does not offer a method for assigning uncertainties to Manning's coefficients, the value ranges provided in the reference tables represent a good starting point. Additional research is needed to verify the PEM analysis results, by comparing the results of flood estimates obtained with a 2D hydraulic model to the results obtained by more accurate methods, such as first-order second-moment (FOSM) and Monte Carlo analysis. Some simplification hypotheses were considered in the analyses. The hypothesis of independence between the input variables was considered to facilitate the estimation of the mean and standard deviation as descriptors of the uncertainty.
The results obtained with the PEM should be considered as approximate values and it should be kept in mind that this method provides a rough estimate of the uncertainty. It is an interesting alternative to computationally more demanding methods in terms of time and money in the risk analysis context. The method is simpler to apply, requiring less information to describe the model inputs and fewer model executions and computations to calculate uncertainty. The method required only the mean and variance of the model inputs in order to complete the analysis.
The uncertainty analysis presented in this study provides a general method for better understanding the H2D2 model. It provides useful information for researchers and managers, particularly if computations are conducted in near real time during an inundation event. It can raise awareness of the main uncertainties of models.