Uncertainty Analysis of Spatiotemporal Models with Point Estimate Methods (PEMs)—The Case of the ANUGA Hydrodynamic Model

Abstract

Practitioners often neglect the uncertainty inherent to models and their inputs. Point Estimate Methods (PEMs) offer an alternative to the common, but computationally demanding, method for assessing model uncertainty, Monte Carlo (MC) simulation. PEMs rerun the model with representative values of the probability distribution of the uncertain variable. The results can estimate the statistical moments of the output distribution. Hong’s method is the specific PEM implemented here for a case study that simulates water runoff using the ANUGA model for an area in Glasgow, UK. Elevation is the source of uncertainty. Three realizations of the Sequential Gaussian Simulation, which produces the random error fields that can be used as inputs for any spatial model, are scaled according to representative values of the distribution and their weights. The output from a MC simulation is used for validation. A comparison of the first two statistical moments indicates that Hong’s method tends to underestimate the first moment and overestimate the second moment. Model efficiency performance measures validate the usefulness of Hong’s method for the approximation of the first two moments, despite the method suffering from outliers. Estimation was less accurate for higher moments but the moment estimates were sufficient to use the Grams-Charlier Expansion to fit a distribution to them. Regarding probabilistic flood-inundation maps, Hong’s method shows very similar probabilities in the same areas as the MC simulation. However, the former requires just three 11-minute simulation runs, rather than the 500 required for the MC simulation. Hong’s method therefore appears attractive for approximating the uncertainty of spatiotemporal models.

1. Introduction

Flood inundation, in its many forms, is one the most devastating types of natural disaster for civilization [1]. Floods are involved in the majority of fatalities associated with natural disasters [2] and cause severe economic damage by disrupting and destroying processes within the society and economy. The socio-economic impact of flood inundations has given rise to studies on the subject receiving widespread interest. Flood-risk assessments are in particular focus and are legally required in some parts of the world, for example in the European Directive on the Assessment and Management of Flood Risks [3].

Maps are one of the best ways to present the spatial distribution of uncertainty and risk. Flood-risk maps are event-based and show the probability of the occurrence and consequences of a flood [1]. The process of producing flood-risk maps involves floodplain mapping and uncertainty analysis. Floodplain maps are generated by a deterministic approach using either empirical, hydrodynamic or simplified-conceptual models [4,5,6,7,8,9]. Historical flood data is employed to validate models and calibrate their parameters [10]. Despite a notable research effort in recent decades, deterministic models remain afflicted by uncertainty [9]. This is in part because the models are assumed to be a comprehensive representation of the relevant physical processes and their relationship to water flow, and in part because they are assumed to use an optimal set of parameters. The issue is further compounded because uncertainty is likely to be non-stationary. Practitioners such as environmental agencies, river basin authorities and engineering consultancies often fail to fully invest the effort required for the detailed setting up of model and lack the computational resources required to solve the physical models [1]. This means that, in practice, additional uncertainty can be introduced to flood inundation modelling by the application of simplified, conceptual models in a deterministic way.

Realistically evaluating the risk of flood events and planning countermeasures therefore requires probabilistic floodplain maps, which are produced from ensemble simulations, like sensitivity analysis, of deterministic models [1]. The probabilistic approach is also more robust concerning the non-stationary aspect of uncertainty because of the many simulations that use a diverse set of data and parameters. Despite awareness of the benefits of this approach, Di Baldassarre et al. [1] claims that practitioners rarely apply probabilistic methods because of a slow rate of knowledge transfer and their familiarity with deterministic results. The authors further conclude there is a lack of a coherent terminology, a systematic approach and mature guidance for the use of probabilistic methods in hydrological analyses.

Reliability engineering has made a major contribution to uncertainty analysis. For one-dimensional models, several approaches can be used to quantify uncertainty. These include the First-Order Second-Moment (FOSM) method [11], Monte Carlo (MC) simulations [12], the Moment method [13], the spectral method [14], the Mellin transformation [15,16] and Point Estimate Methods (PEMs) [17,18,19,20,21,22].

Two-dimensional models are complicated by the need to consider spatial variability, which is handled by random field theories. Random fields can be generated in many ways [23] and serve as input for MC simulations which are common in uncertainty analysis, though very demanding in terms of computational resources. The study by Aerts et al. [24] is a good example: the researchers attempt to optimally locate and design a ski run while considering the uncertainty arising from an uncertain Digital Elevation Map (DEM). The uncertainty is analyzed by a MC simulation, the distinct DEM realizations are derivatives of random error fields, and the error values are the result of measurements at ground control locations. The measured differences are then analyzed for their spatial structure in the form of a variogram. Kriging is used as part of a Sequential Gaussian Simulation (SGS) [25,26,27]. SGS is a method of generating random fields to interpolate the non-measured errors between the ground control locations.

The popularity of MC simulations is illustrated by their continuous development. This development has led to various stratified sampling methods [28,29] that reduce the computational burden of MC, for example, Latin Hypercube Sampling [30,31,32]. Another improvement was the combination of MC with the Markov Chain process [33,34,35]. A useful implementation of MC can be found in the work by Feinberg and Langtangen [36].

Suchomel et al. [37] present a hybrid method combining the FOSM with random fields. Although the method can only consider uncorrelated uncertain variables, spatial correlation is accounted for by spatial averaging, which thus implies a variance reduction. The method was implemented as a finite element model.

In hydrology, uncertainty is commonly quantified by sensitivity analysis [9]. This involves perturbing within feasible limits variables and parameters–either individually and gradually (local) or simultaneously and randomly–and then assessing the effects on the outputs and any interactions.

The Generalized Likelihood Uncertainty Estimation (GLUE) is based on a Bayesian approach [1]. Similar to MC, GLUE requires several model realizations using random samples of the parameter in focus. The realizations are subsequently distinguished into behavioral and non-behavioral groups. The outputs (e.g., the extent of flood inundation) are merged using Bayesian equations to result in a distributed probability profile. GLUE is less rigorous and requires the user to make many choices, increasing the potential for methodological errors [9].

A probabilistic modelling approach usually involves repeated runs of deterministic models using samples from across the range of parameters and available data. The purpose is to approximate the probability distribution of the output. For users with scarce computational resources, simplified conceptual methods are often more convenient than physics-based models.

The purpose of this study is to expand the application of PEMs to uncertainty analysis of spatiotemporal models and to develop a convenient way to produce probabilistic inundation maps.

2. Materials and Methods

This section explains the modelling approach and data used in this study. First, the ANUGA inundation spatial model [38,39] is briefly explained. The MC approach is subsequently discussed and this is followed by an introduction to Hong’s PEM [20]. Furthermore, the expansion of PEM to two-dimensional space is clarified. Finally, the section concludes by describing a case study and its data set.

2.1. ANUGA

ANUGA is developed by the Australian National University (ANU) and Geoscience Australia (GA). The fluid dynamics of ANUGA is implemented using a finite volume approach to solve the Shallow Water Wave Equations. The simulation domain is partitioned into a mesh of triangular cells. The governing Shallow Water Wave Equations are solved at the centroid of the cell and thus provide water depth and horizontal momentum over time. ANUGA is capable of simulating the wetting and drying process when water enters and leaves a cell. It can therefore approximate the flow onto a beach, dry land and around structures such as buildings. Discontinuities, such as hydraulic jumps, are accommodated by the finite-volume implementation. A scenario set-up requires the domain’s geometry, initial water level and boundary conditions. The drivers of the system in the inundation model, like rainfall, abstraction and culverts, are designed as operators and are also required to set up the scenario. The mesh generator can be used interactively to help the user to tailor the domain’s geometry. Mesh triangles are symbolically tagged to indicate boundary regions or regions with different parameters, such as the Manning friction coefficients. ANUGA is mainly written in the Python programming language, which allows for flexible usage [39]. Computationally intensive components are executed as C routines by the means of Python’s numerical module Numpy. The approach has been validated by adequately modelling various inundation problems [40,41].

2.2. Monte Carlo Simulation

MC simulation involves running a simulation many times with different random samples from the probability distribution of an uncertain input variable to obtain a numerical probability distribution of the output variable [24]. MC is very commonly used in uncertainty quantification to capture a system’s natural variability [12,42,43]. The MC simulation serves as benchmark in this study to test the accuracy of the spatially implemented PEM. The outcome of MC is currently the best approximation to the “ground truth” of a stochastic system.

Distinct spatial realizations of the input variable in question (ground elevation) are produced by a SGS. Random error fields induce noise into the DEM data to derive distinct realizations. The error values originate from the differences between the original data set and the ground control data set at 50 locations. Figure 1 shows the histogram of the elevation error, which displays a non-normal form. Error values were interpolated by Ordinary Kriging to obtain error values for the whole spatial domain. Because Kriging provides an estimate of a variable’s mean and standard deviation, the variable can be represented as a random variable according to a Gaussian distribution. However, Figure 1 shows that the errors are not normally distributed, and thus need to be normalized by the means of normal score transformation [44].

The MC approach requires many simulation runs. Accordingly, 500 were carried out for this study, requiring 500 distinct realizations of the random error field. The first step was to create the same number of random paths visiting every grid cell. Following these paths, the error value for each non-ground control cell was interpolated according to neighboring cells in relationship to the variogram. Because each path is different, the error value of each cell differed slightly from those around it. The mean and standard deviation of the Kriging enabled a random selection from a normal distribution. SGS was conducted by a Python library using the mean of the High Performance Geostatistics Library [45]. Subsequently, the error values were back-transformed into the original distribution [44]. The advantage of using Kriging was weighting the influence of the neighboring cells according to the semivariogram permitted the spatial structure to be respected.

2.3. Hong’s Point Estimate Method

2.3.1. Theory

The core idea of the PEM is to calculate the statistical moments of the output variables. The moments are derived from several realizations of a deterministic model using representative values of the uncertain variable as an input. The purpose is to approximate the probability distribution of the model’s output variable, which is related to the probability distribution of the uncertain input variables. The process only requires a few statistical moments of the input variables to be known. The derivation of PEM assumes that the rth moment of the output variable is the integral of the evaluation function with the random input variables over their domains with respect to their joint distribution functions. The equations of the representative values and weights are the products of expanding this function by Taylor series at various points of the variable distributions. The original PEM was proposed by Rosenblueth [21], but this has since proved to be too inefficient because it requires $2^m$ model evaluations, where m is the number of uncertain variables. Thus, the refinement suggested by Hong [20], which is a is a $2m+1$ PEM scheme, was used.

Although this study considers only one variable (elevation), the set of random variables is defined as $X_i,i=1,2,\cdots ,m$. Therefore, to compute the statistical moments of the output requires just three simulation runs, one for each representative input value, k. The first four statistical central moments (mean, variance, skewness and kurtosis) were used to approximate the distribution of the uncertain variable and thus compute the location of the representative values.

The locations of the representative values for each uncertain variable, i, in standardized space, ${\xi }_{i,k}$, were computed by Equation (1) based on the skewness, ${\lambda }_{i,3}$, and kurtosis, ${\lambda }_{i,4}$, of the uncertain variable.

$${\xi }_{i,k}=\left\{\begin{array}{cc}\frac{{\lambda }_{i,3}}{2}+{(-1)}^{3-k}\sqrt{{\lambda }_{i,4}-\frac{3}{4}{\lambda }_{i,3}^2}\hfill & ,k=1,2\hfill \\ 0\hfill & ,k=3\hfill \end{array}\right.$$

(1)

Subsequently, the standardized location was transformed using the mean, ${\mu }_i$, and standard deviation, ${\sigma }_i$, of the uncertain variable to its original space, $x_{i,k}$, by Equation (2). The representative value $k=3$ was the mean value of the variable.

$$x_{i,k}={\mu }_i+{\sigma }_i{\xi }_{i,k}$$

(2)

Based on the standardized values for location, skewness and kurtosis, the weights for the representative values, ${\omega }_{i,k}$, were then calculated by Equation (3).

$${\omega }_{i,k}=\left\{\begin{array}{cc}{(-1)}^{3-k}\frac{1}{{\xi }_{i,k}({\xi }_{i,1}-{\xi }_{i,2})}\hfill & ,k=1,2\hfill \\ \frac{1}{m}-\frac{1}{{\xi }_{i,4}-{\xi }_{i,3}^2}\hfill & ,k=3\hfill \end{array}\right.$$

(3)

The weight determines the proportion of impact that each evaluation of the model $F(.)$ at the location of the representative value has on the statistical moment of the output variable. The rth statistical moment of a function of m random variables, Z, was then approximated by Equation (4).

$$\begin{array}{cc}\hfill E\left(Z^r\right)\approx & {\Sigma }_{i=1}^m{\omega }_{i,3}{\left(F({\mu }_{x1},{\mu }_{x2},\dots ,{\mu }_{xi},\dots ,{\mu }_{xm})\right)}^r+\hfill \\ \hfill & {\Sigma }_{i=1}^m{\Sigma }_{k=1}^2{\omega }_{i,k}{\left(F({\mu }_{x1},{\mu }_{x2},\dots ,x_{i,k},\dots ,{\mu }_{xm})\right)}^r\hfill \end{array}$$

(4)

2.3.2. Spatial Implementation of PEM

Similar to the MC approach, the spatial implementation also relies on random error fields. Instead of 500 error fields, this study only uses $2m+1$ error fields, i.e., three, given that only one uncertain variable is considered. The three error fields were produced by the SGS process. Hong’s method was then used to determine the representative values of the distribution as described by Equations (1) and (2). Figure 1 shows the histogram of the errors with the locations of the representative values.

The representative values and the corresponding weights were thereafter used to scale the three error fields, as can be seen in Equation (5).

$$sRE{F}_{i,k}=x_{i,k}-(RE{F}_i-{\mu }_{RE{F}_i})\frac{{\omega }_{i,k}{\sigma }_i}{{\sigma }_{RE{F}_i}}$$

(5)

where $RE{F}_i$ is the SGS-produced random error field of the uncertain variable, i. Figure 2 shows the realizations of the random error field, $sRE{F}_{i,k}$, scaled by the representative values, $x_{i,k}$; their weights, ${\omega }_{i,k}$; and the standard deviation of the uncertain variable, ${\sigma }_i$. It is important that the range of values in three scaled random error fields is similar to that of the original random error field.

The simulation was run three times with the corresponding error fields. The statistical moments of the output variable were computed with Equation (4) as the weighted sum of the simulation results at the locations of the representative values to the power of the order of the statistical moment.

2.4. Workflow

The Figure 3 illustrates the workflow of this study. The point of departure is the DEM and elevation of ground control locations. Both data sets combined result in error values. The error values need to be normalized because SGS requires that the input variable has a normal distribution. SGS generates 500 error fields for MC and three for PEM. The error values are back-transformed to the original distribution. Equation (5) scales the error fields for PEM. The error fields modify the DEM. The modified DEMs are supplied to ANUGA with the purpose to run an inundation scenario with a rainfall operator. The outputs are 500 raster time series data sets for MC and three for PEM. In the case of MC, the 500 raster time series data sets are reduced to four data sets, one for each statistical moment. In the case of PEM, Equation (4) estimates the raster time series data of the moments. The model efficiency performance measures finally compare the raster time series of each moment estimated by both methods. For each model efficiency performance measure, the results are one raster data set for each moment where each grid cell contains a model efficiency value. Each raster data set is subsequently summarized by mean, median, standard deviation, range and skew. The flowchart only displays the model efficiency performance assessment for the reason of readability.

2.5. Case Study

The location of the case study is Cockenzie Street and its surroundings in the city of Glasgow, UK. The data set was provided by the UK Environment Agency following their benchmark test No. 8 of two-dimensional hydraulic packages [41]. The data set contains the DEM, the geometry of buildings and roads, the Manning friction coefficients of paved and unpaved surfaces, the extent of the simulation domain, the location of the nine water level gauges and a rainfall curve. The study area is visualized in Figure 4. The DEM covers an area of approximately 0.4 km by 0.96 km. The original resolution of the DEM (0.5 m) was reduced to 5 m because of computational limitations. The range of the ground elevation is 21 m to 37 m. There are over 14.000 grid cells in the domain.

There are only two land cover categories with Manning friction coefficients: pavements (0.02) and unpaved surfaces (0.05). The simulation time was 1000 s with a 10-second time-step. The initial condition was a dry bed. The rainfall began after 1 min and lasted for 3 min with an intensity of 400 mm/h.

3. Results and Discussion

The model comparison below follows the recommendations in Bennett et al. [46] and Biondi et al. [47]. The times series at the nine gauge locations indicated in Figure 4 were visually compared. The similarities and differences of the spatial extents of the maximal inundations are then briefly assessed. The evaluation is completed by using the following model efficiency performance measures [47] on all time series over the whole spatial domain: Nash-Sutcliffe Model Efficiency (NSE), Root Mean Square Error (RMSE) and Index of Agreement (IoA). The statistical moments of MC should correspondingly be considered as “ground truth”, i.e., akin to observed data. We focus on the first and second statistical moments because of their more important relationship with real-world applications, but also consider the third and fourth moments because they permit the fitting of a probability distribution to the output variable.

A comparison of the water depth values is shown using residual plots in Figure 5 and Figure 6. The residual plots show the differences and similarities over time. Figure 6 groups the time series with noteworthy behavior. During rainfall (from time-step 61 to 240), PEM mostly underestimates the first and second moment, except for at gauges 4 and 6 where it overestimates both moments. Gauge 7 is particular as PEM underestimates the first moment here but then overestimates the second moment. PEM also appears to be less sensitive during the more dynamic period, i.e., when mass (rainfall) is added to the system. Therefore, the residuals are much higher during these periods. After the rainfall had stopped, the output from MC and PEM converged at most gauges, except at gauges 1 and 3. PEM continued to slightly underestimate the first moment at most gauge locations while also overestimating the second moment at many gauge locations. Concerning the first moment, there is a bias towards MC, except at gauge 4. The plots show a bias of the second moment towards PEM, which is especially visible for gauge 3. The estimations of the second moment by PEM are also marked by outliers, especially for gauges 2, 6, 7 and 9. However, it is very sensitive to the multitude of paths runoff can take when slightly changing elevation. This is especially the case for the gauges grouped in Figure 6.

The spatial extent of the first two moments of the maximal flood inundation are compared in both Figure 7 and Figure 8 because peak water depth is important for practical purposes. The mean estimations by MC and PEM in Figure 7 depict a very similar spatial structure and magnitude. Figure 8 confirms the previous statements about the second moment. Despite similar spatial structures, the second moment is systematically overestimated by PEM. This observation is useful as it seems that PEM sets the upper limit for the second moment.

Table 1. Summary statistics of model efficiency performance for all time series of the estimated statistical moments at each grid cell.

	${\mathit{NSE}}_{\mathit{mean}}$	${\mathit{NSE}}_{\mathit{median}}$	${\mathit{NSE}}_{\mathit{SD}}$	${\mathit{NSE}}_{\mathit{range}}$	${\mathit{NSE}}_{\mathit{skew}}$
1st Moment	0.73	0.99	2.35	141.49	−34.87
2nd Moment	−1.95	0.96	38.10	2008.66	−35.86
3rd Moment	−0.78	−0.25	1.81	12.61	−1.72
4th Moment	0.07	−0.32	0.96	10.92	−1.74
	${\mathit{RMSE}}_{\mathit{mean}}$	${\mathit{RMSE}}_{\mathit{median}}$	${\mathit{RMSE}}_{\mathit{SD}}$	${\mathit{RMSE}}_{\mathit{range}}$	${\mathit{RMSE}}_{\mathit{skew}}$
1st Moment	0.0074	0.0007	0.0164	0.1775	3.6896
2nd Moment	0.0179	0.0014	0.054	0.9263	6.1075
3rd Moment	5.092	3.6726	4.101	20.4579	1.0632
4th Moment	74.6667	33.7397	85.76	443.706	1.7032
	${\mathit{IoA}}_{\mathit{mean}}$	${\mathit{IoA}}_{\mathit{median}}$	${\mathit{IoA}}_{\mathit{SD}}$	${\mathit{IoA}}_{\mathit{range}}$	${\mathit{IoA}}_{\mathit{skew}}$
1st Moment	0.96	0.99	0.09	0.96	−5.28
2nd Moment	0.88	0.99	0.18	0.99	−2.22
3rd Moment	0.55	0.64	0.35	0.99	−0.23
4th Moment	0.62	0.81	0.36	0.99	−0.51

summarizes the model efficiency performance measures for all time series of the estimated statistical moments. The model efficiency performance measures consequently compare the similarity between the results of MC and PEM. The empirical distributions of the performance measures are described by mean, median, standard deviation, range and skew. The column $NS{E}_{mean}$ characterizes the means of the Nash–Sutcliffe efficiency between MC and PEM of every estimated moment. The other columns follow the same logic. The NSE in Table 1 is the most commonly employed performance measure in hydrology. It ranges from 1 to negative infinity, where 1 indicates a perfect match. The RMSE in Table 1 ranges from 0 to positive infinity, where 0 indicates a perfect match. The NSE and RMSE suffer potential bias when comparing models with offsets [46], perhaps because of the usage of the square root. To complement the evaluation, the IoA in Table 1 was also calculated. The advantage of using the IoA is that it uses absolute differences, and thus might better avoid bias. It ranges from 0 to 1, where 1 indicates a perfect match. The mean values for the performance measures for the first and second moments show a good match. The situation is further improved when regarding the median owing to the impact on the mean values of far outliers, as indicated by the range and skew values. The outliers also prevent the performance measures being summarized in a box plot. Estimations of the second moment by PEM yielded more under-outliers, as is clearly observed when comparing the negative mean of the NSE with the median. The third and fourth moments perform less well. The reasons for this appear to be more systematic because range and skew indicate that outliers are less of a problem than for the first and second moments. Interestingly, the NSE and IoA suggest that the estimations of the fourth moment match the observations better than those of the third moment, but RMSE suggests the opposite. The decreasing performance of the PEM estimations with the order of moments might have been caused by the usage of the square function in the calculation of the moments in Equation (4).

The end-product of this study is the probabilistic inundation map. The Gram-Charlier expansion, which was implemented using the R package PDQutils [48], was used to fit the statistical moments to a probability distribution based on a normal distribution. Figure 9 displays the probability of a maximal inundation over 0.2 m. Despite the differences of the higher moments between MC and PEM, PEM predicts a similar risk of inundation of over 0.2 m in terms of probability values and spatial extent.

4. Conclusions

This study demonstrated that PEMs are capable of approximating the uncertainty of an output in spatiotemporal models. Using PEMs instead of the more commonly used MC simulations yields considerable savings of computational resources. Instead of 500 simulation runs, each lasting 11 min for the data used in this study, only 3 simulation runs were required for the PEM approach. This decrease in computational resource requirement is especially important for spatiotemporal models because of the amount of data and the complexity involved. Methodological advances like this are key to development because challenges arising from the use of Big Data are unlikely to be overcome by increases in computational power alone, despite its continually falling cost.

Uncertainty analysis is important to mitigate the impact of an overreliance by practitioners on deceptive deterministic solutions that arise because fundamental natural principles are less deterministic than we think and our continued lack of capacity to model phenomena in their totality. While MC simulations are computationally too demanding to be widely adopted outside of the research community, PEMs’ ability to decrease computational costs could facilitate the uptake of uncertainty analysis in various fields. For example, hydrologists could combine water pollution models and PEMs to assess the risk of exceeding the legislated values (such as the U.S. EPA’s Total Maximum Daily Loads). In a similar fashion, PEMs could help to determine the risk of exceeding health-related air pollution thresholds. Concerning flooding, flood-control authorities and insurance companies could better evaluate the potential risks for critical infrastructure and for insured objects, respectively.

Future studies should investigate different implementations of PEMs for spatiotemporal models and how to accommodate multivariate cases, where spatial cross-correlation is of tremendous importance. Current visualization techniques need to be improved to exploit the full potential of the information furnished by uncertainty analysis.