Uncertainty Quantification and Sensitivity Analysis Applied to the Wind Wave Model Swan

(1)

UNCERTAINTY QUANTIFICATION AND SENSITIVITY ANALYSIS APPLIED TO THE WIND WAVE MODEL SWAN

By

Anna Nikishova

Thesis submitted to

Graduate School of Informatics Computational Science Lab

University of Amsterdam and

Faculty of Informational Technologies and Programming Department of High Performance Computing

Saint Petersburg National Research University of Information Technologies, Mechanics and Optics

In partial fulfillment of the requirements for the degree of

Master of Computational Science

Supervisors: Prof. Alfons Hoekstra Dr. Valeria Krzhizhanovskaya

Scientific advisers: Dr. Sergey Kosukhin Prof. Alexander Boukhanovsky

Amsterdam, Saint Petersburg 2016

(2)

2

STATEMENT OF ORIGINALITY

This document is written by Student Nikishova Anna who declares to take full responsibility for the contents of this document.

I declare that the text and the work presented in this document is original and that no sources other than those mentioned in the text and its references have been used in creating it.

(3)

3

The thesis is dedicated

to my mother Elena and Giovanni Eugenio Comi

for your constant support and love,

(4)

4

ACKNOWLEDGEMENTS

My sincere gratitude to my supervisors Prof. Alfons Georgius Hoekstra and Dr. Valeria Vladimirovna Krzhizhanovskaya for guiding me in this research.

I would like to thank my scientific advisors Dr. Sergey Sergeyevich Kosukhin and Prof. Alexander Valerievich Boukhanovsky for inspiring me to carry out this work and their help in the experiments design in this project.

I am deeply grateful to Dr. Giovanni Eugenio Comi for his support and helping me with some mathematical aspects lying on the basis of this research.

I would like to extend my appreciation to my teachers for sharing with me their knowledge and skills, without which I would not be able to perform this research.

(5)

5

ABSTRACT

Keywords: uncertainty quantification, sensitivity analysis, Sobol method, Monte Carlo simulations, wind wave model, SWAN, Baltic Sea, unstructured mesh.

The purpose of this study is to quantify uncertainty in the response from the computations of the wind wave model SWAN on unstructured mesh, and to investigate the effects of uncertainty in the input parameters to the overall model uncertainty.

We applied the Bayesian network approach in order to estimate model uncertainty. The Sobol method was used to analyze input parameters contribution to uncertainty in the model response. Both the methods are based on the Monte Carlo approach. We used the low-discrepancy Sobol sequence in order to decrease the number of samples and to obtain more accurate results.

In the thesis, we show that the most important input parameter is wind velocity. In addition, sensitivity to uncertainty in bottom level values depends on location. Since unstructured meshes allow represent spatial data with different resolution in different areas, we are able to increase mesh resolution in more sensitive regions, and, in such a way, decrease uncertainty in the model response. Some more details about regular and unstructured computational grids for surge flood forecasting and prevention we presented in (Kosukhin, et al., 2015).

(6)

6

INTRODUCTION

Applications of computational modelling become more and more popular in science and technology. In some cases, the cost of computer simulations is incommensurably less than a physical experiment. Alternatively, it can be simply not possible or, for example, not ethical to execute an experiment and, in such cases, a computer model substitutes physical implementation of the process.

The goal of computational experiments is usually either to understand and gain knowledge of a process, the so-called “exploratory” modelling, or to predict the behaviour of the system, i.e. the “simulation” modelling. Important stages of any computer simulation are model validation and verification (Oberkampf, et al., 2002).

To verify a model we need to evaluate the accuracy through matching the output from simulation with solutions known to be highly accurate (Figure 0.1). The goal of the model verification is to check whether the model is built in a right way (Balci, 1995). There exist different techniques to verify models, however, we did not apply them in this work, since we were concerned with the study of model validation, whose definition is given below.

Figure 0.1 – The scheme of verification processes (AIAA (American Institute of Aeronautics and Astronautics), 1998)

(8)

8

The model validation investigates if the model produces the result with required accuracy. During the validation process, it is analyzed whether the right model was built (Balci, 1995). The output from computational model is compared and analyzed together with experimental data from various sources (Figure 0.2). Fundamental components of model validation are uncertainty quantification and sensitivity analysis.

Figure 0.2 – The scheme of validation processes (AIAA (American Institute of Aeronautics and Astronautics), 1998)

Uncertainty quantification is the science of estimation the distribution of the result values from computational simulation in conditions of a lack of knowledge about exact value of input parameters. Uncertainty quantification can be used to evaluate and characterize the behaviour of the model computations and the likelihood of different outcomes (Crosetto, et al., 2000).

Sensitivity analysis allows one to estimate the influence of uncertain input on the variability of model output. The increase of interest to Sensitivity Analysis (Figure 0.3) was demonstrated in (Song, et al., 2015).

(9)

9

Figure 0.3 – Number of publications on Sensitivity Analysis in corresponding year. The dashed bars represent the search in all articles; the white bars denote the result based on papers only from ‘‘water resources’’ and ‘‘environmental sciences’’ categories (Song, et

al., 2015)

The aim of this study was to quantify the uncertainty and apply sensitivity analysis to the simulation model SWAN (Simulating WAves Nearshore) for wind wave predictions for the Baltic Sea region. The SWAN model (SWAN Team, 1993-2015) is a third-generation numerical wave model to compute random, short-crested waves in coastal regions with shallow water and ambient currents. The observed input parameters were bottom level, bottom friction coefficient, wind velocity, whitecapping dissipation coefficient and wave steepness coefficient. The output parameter under consideration was significant wave height. We used an unstructured mesh, defined over the Baltic Sea with increased resolution in the area of the Gulf of Finland. In Chapter 1, we discuss the model in some more details, its configuration and the unstructured mesh used in this study.

(10)

10

Chapter 2 is devoted to the notions of uncertainty in computational modelling and uncertainty quantification using the Bayesian network approach. This methodology allows establishing statistical relationships between input factors and model output using the principle of uncertainty propagation through the network with the aid of Monte Carlo methods. The advantage of the method is that there is no need to perform any modifications of the model, and the approach can be applied to computational models with a complex hierarchical structure (Sankararaman, et al., 2012).

The Sobol method (Sobol, et al., 2005) was applied to establish the effect of uncertainty in an input factor on uncertainty in model response. The method is based on the decomposition of variance, introduced by Ilya Meyerovich Sobol in 1990. Variance based approaches are effective and flexible methods, because such methodologies do not deal with solvers inside the model and can be applied to models with unknown or complex structure. The detailed representation of the Sobol method is given in Chapter 3. In the thesis, quasi-Monte Carlo method (Sobol, 1998) was implemented using low-discrepancy Sobol sequence. The description and advantages of this approach are shown in Chapter 4.

The results of uncertainty quantification and sensitivity analysis are provided in Chapter 5. The main investigation of the research was the detection that the wind parameter is the most influential among ones included in the study. Conclusions and future work are presented in the last chapter of the thesis.

(11)

11

1 Wind wave model SWAN

The city of Saint Petersburg is threatened by floods, because large parts of the city are slightly above sea level. Starting from the time when Saint Petersburg was founded in 1703, flood conditions have occurred more than 300 times. Since 2011, Saint Petersburg is protected by the Flood Prevention Facility Complex. The gates of the prevention complex (Figure 1.1) can be closed in case of a warning situation in order to hold back water. The decision about gates operating are based on forecasts, obtained by coupling models and data from different sources. In this thesis, we applied uncertainty quantification and sensitivity analysis methods to the wind wave model SWAN, which is a part of the models ensemble.

Figure 1.1 – The Saint Petersburg Flood Prevention Facility Complex (Kosukhin, et al., 2015)

The third-generation numerical wave model SWAN (Simulating WAves Nearshore) has been developed at Delft University of Technology, the Netherlands. The model simulates random, short-crested waves, in order to make predictions about wave parameters in coastal areas, lakes and estuaries (Booij, et al., 1999).

(12)

12

Five model parameters were under consideration: wind velocity, bottom level, bottom friction coefficient, whitecapping dissipation coefficient, and wave steepness coefficient. The representation of this factors in the model is shown below. A more detailed description of model parameters can be found in (SWAN Team, 1993-2015).

1.1 Model input parameters

The model inputs are wind velocity and bottom level map. Traditionally, the spectral wave energy equation is the basis of wave models (Tolman, 2008). In SWAN, all physical processes of wave energy are described by the following equation:

𝑆_𝑡𝑜𝑡 = 𝑆_𝑖𝑛+ 𝑆_𝑛𝑙3+ 𝑆_𝑛𝑙4+ 𝑆_{𝑑𝑠,𝑤} + 𝑆_{𝑑𝑠,𝑏} + 𝑆_{𝑑𝑠,𝑏𝑟}, (1.1) where, 𝑆_𝑖𝑛 denotes wave growth by wind, 𝑆_𝑛𝑙3 and 𝑆_𝑛𝑙4 are nonlinear wave energy transfer from three-wave and four-wave interactions, 𝑆_{𝑑𝑠,𝑤}, 𝑆_{𝑑𝑠,𝑏} and 𝑆_{𝑑𝑠,𝑏𝑟} are wave decay due to whitecapping, bottom friction and depth-induced wave breaking, respectively. The right hand side of equation (1.1) represents sources and sinks of wave energy.

1.1.1 Wind velocity

In SWAN, wind speed is defined at 10 meters above sea level, and the value should be given by two components, which correspond to a two-dimensional horizontal coordinate system, defined by the user. In (1.1), The wind wave growth is given by

𝑆_𝑖𝑛(𝜎, 𝜃) = 𝐴 + 𝐵𝐸(𝜎, 𝜃), (1.2) where 𝜎 and 𝜃 are wave frequency and direction, respectively, E is wave energy, 𝐴 denotes linear growth, and 𝐵 describes exponential growth.

The values of wind speed were taken using simulations with the forecasting model HIRLAM – HIgh Resolution Limited Area Model (HIRLAM, 2002). For nonstationary runs of SWAN, wind data were taken with a time step of six hours.

(13)

13

1.1.2 Bottom level and bottom friction coefficient

The bottom level map contains information about the water depth for the area of interests. In this study, the bathymetry data was given on an unstructured computational grid presented in the part 1.3.

For sandy bottom continental shelf seas, bottom friction is the dominant mechanism, which has the form of

𝑆_{𝑑𝑠,𝑏}(𝜎, 𝜃) = −𝐶𝑏

𝜎2

𝑔2_sinh2_𝑘𝑑𝐸(𝜎, 𝜃), (1.3)

where 𝐶_𝑏 is the bottom friction coefficient, 𝑔 is the gravity acceleration, 𝑘 is wave number, and 𝑑 is depth. In the model, the friction coefficient can be a constant value or represented on a spatial grid. We defined the value as a constant using the JONSWAP option (Hasselmann, et al., 1973).

1.1.3 Whitecapping dissipation rate and wave steepness coefficients

The whitecapping processes is described by the Hasselmann pulse-based model (Hasselmann, 1974):

𝑆_{𝑑𝑠,𝑤}(𝜎, 𝜃) = −Γ𝜎̃𝑘

𝑘̃ 𝐸(𝜎, 𝜃), (1.4) where value of Γ dependents on wave steepness coefficient and coefficient for defining the whitecapping dissipation rate, 𝑘̃ and 𝜎̃ denote mean wave number and mean frequency, respectively.

1.2 Model output parameters

The model produces such wave parameters as significant wave height, mean wavelength, mean absolute wave period, peak direction, direction of energy transport and

(14)

14

many other (see (SWAN Team, 2013) for more information). The model uncertainty and sensitivity were evaluated with respect to significant wave height 𝐻_𝑠 defined as

𝐻_𝑠 = √4𝐸, (1.5) 𝐸 = ∫ ∫ 𝐹(𝜎, 𝜃) 𝑑𝜎 𝑑𝜃 , (1.6) where F is the wave energy spectrum.

1.3 Model configurations

In order to increase the accuracy of the forecast, the walls of the Flood Prevention Complex were represented in the model using the bathymetric feature “obstacle”. Figure 1.2 illustrates how the walls were denoted in the model.

Figure 1.2 – Green lines indicate the location of the obstacles defined in the SWAN

The command “OBSTACLE” contains the option to define a reflection coefficient. Usually, in order to set the value correctly, knowledge about the material of the physical obstacle is required. We have calibrated this parameter using experimental data to estimate the mean square error given by

(15)

15 𝑀𝑆𝐸 = 1 𝑁∑(𝐸𝑖 − 𝑀𝑖) 2_, 𝑁 𝑖=1 (1.7)

where 𝑁 is the number of forecast points, 𝐸_𝑖 and 𝑀_𝑖 are values from experiment and model, respectively, for i-th forecast point. We tested the values of the reflection coefficient in the range from 0.0 to 1.0 with step 0.1. The result mean square error values are plotted in Figure 1.3. As we can see, the error is lower for the values of the reflection coefficient greater than 0.6. As the result, the value of this coefficient was set equal to 0.7.

Figure 1.3 – Mean square error in the value of significant wave height from runs in non-stationary mode with different values of reflection coefficient

1.4 Unstructured mesh

The choice of geographical computational grid plays a crucial role in model simulations. A usual approach is a regular mesh with additional runs on a nested mesh in order to decrease the overall computational demand. However, there is an important issue

(16)

16

that such orthogonal structured meshes are not sufficiently flexible to fit an arbitrary geometry. Results presented in (Zijlema, 2010) and (Mavriplis, 2007) show that computations on unstructured meshes perform stable and more accurate model outcomes for areas of interest. Moreover, with such approach it is possible to decrease computational time by reducing the density of points in remote areas.

The computational grid used in this research is shown in Figure 1.4. Since the area of our interest was the vicinity of the Saint Petersburg Flood Prevention Facility Complex, the resolution of the mesh was increased in this region. In addition, some coastal regions have more dense mesh due their complex coastal topology. The mesh was constructed using the two-dimensional quality mesh generator Triangle (Shewchuk, 2002).

Figure 1.4 – The unstructured computational grid we used in this research

In the next chapters, we present methods for Uncertainty Quantification (Chapter 2) and Sensitivity analysis (Chapter 3) we applied in the work. The results of the study we show in Chapter 5.

(17)

17

2 Uncertainty Quantification

The aim of the validation process is to evaluate the degree to which the model implementation corresponds to some process in nature (Oberkampf, et al., 2002). There exist different sources of uncertainty (Figure 2.1), which lead to uncertain computational outcomes. For instance, uncertainty can come from initial and boundary conditions,spatial and temporal discretization, model structure or numerical inaccuracy, etc.

Figure 2.1 – Sources of uncertainty and propagation processes (Song, et al., 2011)

The sources and propagation of uncertainty in wave simulations are presented in Figure 2.2. In (Slingo, et al., 2011), it is noted that in weather forecasting it is important to separate model uncertainty and uncertainty that comes from different sub-gridscale phenomena, whose representation is not adequate because of the model resolution (Figure 2.3). Some more details about multiscale sources of uncertainty in hydrological modelling are given in (Song, et al., 2011).

Figure 2.2 – Sources of uncertainties in wave model simulations (Ambühl, et al., 2015)

(18)

18

The inherent variation of the physical system under consideration is described by aleatoric uncertainty (from the Latin word alea, meaning "dice"). This kind of uncertainty is usually treated as stochastic and irreducible.

Epistemic uncertainty (from Greek epistēmē, meaning "knowledge, understanding") arises from incomplete information or lack of knowledge. This kind of uncertainty is also known as reducible, i.e. it is possible to decrease epistemic uncertainty by gathering more precise information or by refining the model (Kiureghian, et al., 2009).

(a) (b)

Figure 2.3 – Portioning of uncertainty and how the sources evolve with time in the example of England rainfall for the (a) 2020s and (b) 2080s from the UK Climate

Projections (Slingo, et al., 2011)

Uncertainty Quantification (UQ) is the science of estimation, analysis and reduction of uncertainty in the model response (Karniadakis, et al., 2006). In (Booker, et al., 2011), uncertainty denotes ‘‘what is not known precisely’’, and quantification stand for the process of estimation and analysis based upon data, models, examination, and so on. Figure 2.4 is a

(19)

19

schematic representation of the process of Uncertainty Quantification and Sensitivity Analysis in hydrological modelling.

Figure 2.4 – Uncertainty and sensitivity analysis in hydrological modelling (Song, et al., 2015)

2.1 Methods for Uncertainty Quantification

A wide range of different methods can be applied to the UQ problem. The methods for analyzing and propagating epistemic uncertainty, such as interval analysis, the Dempster-Shafer evidence theory, and second-order probability, are discussed in (Swiler, et al., 2009). A least-squares approach for linear uncertainty analysis is presented in (Cho, et al., 2015). A reliability-based, a classical hypothesis testing, a Bayesian hypothesis testing, and an area metric-based methods are explored in (Ling, et al., 2013).

In general, we can classify the UQ problem into two types: forward uncertainty propagation, where the various sources of uncertainty are propagated through multiple levels of model to predict the total uncertainty in the result of computations, and inverse UQ, where the model’s parameters of multiple subsystems are calibrated simultaneously using test data.

(20)

20

Methods based on uncertainty propagation can be classified into intrusive or non-intrusive, statistical or deterministic. In intrusive methods, a stochastic model substitutes the original one. In contrary, non-intrusive methods do not require any modification of the model, and usually uncertainty is estimated with respect to multiple model computations.

Statistical quantities of interests are estimated in statistical methods. Examples of such method are simple random or more advanced sampling techniques. The Monte Carlo method with quasi-random sampling for UQ is discussed in (Crosetto, et al., 2000). On the contrary, in deterministic methods, model uncertainty is evaluated with respect to local partial derivatives. Examples of the deterministic approaches can be the Propagation of Errors method or the Adjoint Sensitivity Analysis Procedure (Briggs, 2008).

In (Ambühl, et al., 2015), to validate results from wave state simulations, measurements from a buoy at a certain location and the same time period are used: total uncertainty in the computations from the model SWAN and some other popular spectral wave models is evaluated according to the statistical measures like bias

𝐵 = 1 𝑛∑(𝑜 − 𝑚)𝑖 𝑛 𝑖=1 , (2.1) root-mean-square error 𝑅𝑀𝑆𝐸 = √1 𝑛∑(𝑜 − 𝑚)𝑖 2 𝑛 𝑖=1 , (2.1)

and scatter index

𝑆𝐼 = 𝑅𝑀𝑆𝐸

(21)

21

where n is the total number of data points, 𝑚_𝑖 is a result from simulations and 𝑜_𝑖 is the corresponding measurements from a buoy.

In the thesis, we have applied the Bayesian network approach for UQ presented below. The motivation for such choice was a relatively easy implementation, applicability to complex models, and the fact that the method is based on Monte Carlo simulations, which is coherent with the choice of the Sobol method we used for Sensitivity Analysis.

2.2 Bayesian methodology for Uncertainty Quantification

A wide variety of natural processes occurs on multiple scales. Some systems are represented by complex models with multiple levels. Bayesian approach allows for full accounting the overall uncertainty in multiple scales models and the uncertainty in a particular model level ( Ferreira, et al., 2007). As an example, in (Urbina, et al., 2012), (Sankararaman, et al., 2012) and (Sankararaman, et al., 2014), the Bayesian Network (BN) method is applied to estimate and propagate uncertainty in a complex hierarchical system.

The Bayesian Network approach is used to solve both the forward (uncertainty propagation) and inverse (parameters calibration) problems. Monte Carlo methods are usually applied in the forward problem, where all sources of uncertainty at lower levels are propagated through the Bayesian network to evaluate system-level prediction uncertainty. Markov Chain Monte Carlo-based techniques are often applied to solve the inverse problem and calibrate variables at any model levels using data from all model levels (Sankararaman, et al., 2014).

A Bayesian Network (Figure 2.5) is a graphical representation of the connections between model uncertain quantities (de Oude, 2010), where nodes are the model variables and edges define dependence based on conditional probabilities. In order to predict overall model uncertainty, input parameters uncertainty is propagated through the Bayes network.

(22)

22

Figure 2.5 – Example of a Bayesian Network

Bayesian methods are data analysis techniques based on the principles of Bayesian inference. To illustrate the Bayesian philosophy, let us consider a collection of events 𝐴_𝑖, with 𝑖 = 1, . . , 𝑛, from the sample space Ω, and let 𝐵 be another event from the sample space Ω. The law of conditional probability states that

𝑃(𝐴𝐵) = 𝑃(𝐴|𝐵)𝑃(𝐵) = 𝑃(𝐵|𝐴)𝑃(𝐴). (2.4) From this equation, the Bayes theorem follows:

𝑃(𝐴_𝑖|𝐵) = 𝑃(𝐵|𝐴𝑖)𝑃(𝐴𝑖)

∑ 𝑃(𝐵|𝐴_𝑗 _𝑗)𝑃(𝐴_𝑗), (2.5) for any 𝑖 = 1, . . , 𝑛, where 𝑃(𝐴_𝑖) are the prior probabilities, 𝑃(𝐴_𝑖|𝐵) are called posterior probabilities, and the term 𝑃(𝐵|𝐴_𝑖) is so-called “likelihood”, denoted the probability of event 𝐵 in condition 𝐴_𝑖 has occurred. Using this theorem, probabilities of a model state can be updated taking under consideration the new available knowledge.

Let us consider a simple model which is represented by the directed acyclic graph 𝐆 shown in Figure 2.5, where 𝑿 = {𝑋₁, 𝑋₂, 𝑋₃} is a set of uncertain input variables, and 𝒀 = {𝑌₁, 𝑌₂} is a set of model outputs. The set of directed edges 𝑬 = {(𝑋₁, 𝑌₁), (𝑋2, 𝑌1), (𝑋2, 𝑌2), (𝑋3, 𝑌2), (𝑌1, 𝑌2)} represents conditional probability between

variables (de Oude, 2010). In this example, the node 𝑌₁ represents a lower level model result, and a higher level model response 𝑌₂ depends on this output value. We have assumed that the graph 𝐆 satisfies the Markov Condition, which states that for each 𝑋_𝑖 ⊆ 𝑮, 𝑋_𝑖 is

(23)

23

independent from 𝐆\(𝐃𝐞𝐬𝐜(𝑋_𝑖) ∪ 𝐏𝐚(𝑋_𝑖)), where 𝐃𝐞𝐬𝐜(𝑋_𝑖) are the descendants of 𝑋_𝑖 and 𝐏𝐚(𝑋_𝑖) are parents of 𝑋_𝑖. In this case, the whole network is described by the joint probability density function (Jensen, 2001):

𝑃(𝐺) = ∏ 𝑃(𝑋_𝑖|𝐏𝐚(𝑋_𝑖))

𝑖

, (2.6)

where 𝑋_𝑖 are the i-th node of the network. The distributions of the model responses are given by

𝑓_𝑌₁ = ∫ 𝑓_𝑌₁(𝑦₁|𝑥₁𝑥₂)𝑓_𝑋₁(𝑥₁)𝑓_𝑋₂(𝑥₂) 𝑑𝑥₁𝑑𝑥₂, (2.7)

𝑓_𝑌₂ = ∫ 𝑓_𝑌₂(𝑦2|𝑦1, 𝑥2, 𝑥3)𝑓𝑌₁(𝑦1)𝑓𝑋₂(𝑥2) 𝑓𝑋₃(𝑥3)𝑑𝑦1𝑑𝑥2𝑑𝑥3. (2.8)

To obtain new posterior distribution of the model parameters, prior knowledge of model parameters couples with measured data using Bayes’ rule (2.5). The posterior distribution of the model parameters conditioned by the measured data is calculated by

𝑃_{𝑝𝑜𝑠𝑡}(𝑥_𝑖|𝑦_{𝑚𝑒𝑎𝑠}) = 𝑃𝑀(𝑦𝑚𝑒𝑎𝑠|𝑥𝑖)𝑃𝑝𝑟𝑖(𝑥𝑖)

𝑃(𝑦_{𝑚𝑒𝑎𝑠}) , (2.9) where 𝑦_{𝑚𝑒𝑎𝑠} is the measured data, 𝑃(𝑦_{𝑚𝑒𝑎𝑠}) is the probability density of measured data and 𝑃_𝑝𝑟𝑖(𝑥_𝑖) is the prior distribution of i-th model parameter.

We applied the BN method in order to solve forward problem and estimate uncertainty in the response of the wind wave model SWAN, discussed in Chapter 1 using the quasi-Monte Carlo method with Sobol sequences (Chapter 4) for uncertainty propagation through the BN.

(24)

24

3 Sensitivity Analysis

Sensitivity analysis (SA) of model response (Saltelli, et al., 1993) is the study of the relationships between input and outcome. More specifically, with SA methods it is possible to estimate the effects of uncertainty in the input parameters on uncertainty in the model result.

SA techniques can be based on analysis of model behaviour, or can use variance-based methods. Song et al. presented an extensive study on the SA methods (Song, et al., 2015) and the summary of the comprehensive research is shown in Table 3.1.

In the thesis, the Sobol variance-based method for global sensitivity analysis has been applied. The main goal of variance-based methods is to estimate the contribution of the variance in an input factor (or group of factors) with respect to the total variance of the output factor. The advantage of this approach is that there is no need of any assumption about the model structure. Moreover, with these methods, it is possible to distinguish between first-order and higher order effects, and to make corresponding analysis and to draw conclusions. The method is discussed in more detail below.

Along with global approaches, there exists local sensitivity analysis methods. To illustrate the difference, let our model be described by the function 𝑢 = 𝑓(𝑥), where x = (x1,…,xn) are the input parameters, which can be seen as coordinates of a point inside a n-dimension space, and u is the model response, which is a scalar. Local sensitivity analysis is based on partial derivatives of given function. In local methods (Castillo, et al., 2006), only one of possible solutions u*_{= f(x}*_{) is under consideration and the sensitivity of the outcome} u*_{is estimated only with respect to x}*_{. On the other hand, in global sensitivity analysis, there}

is no consideration of individual solutions. Global methods, e.g. the Sobol method, are based on the ANOVA representation (or decomposition), where the abbreviation “ANOVA” refers to the analysis of variance.

(25)

25 Table 3.1 – Comparison of the SA methods (Song, et al., 2015)

Methods Description Application

Local

Estimate local sensitivity measures. Consider model response in a single location in the parameter space. The methods are easy in

realization and interpretation, have low computational cost.

Local analysis of input parameters. Global The main and join effects evaluated over the entire space of uncertain

parameters. Require multiple model evaluations. Identification of main and join effects.

Mathematical Investigate local or linear effects according to single parameter. Deterministic analysis of linear models. Verification and validation. Statistical Based on sampling methods with multiple model evaluations. Probabilistic analysis. Verification. Graphical The methods are used in addition to the mathematical or statistical

methods to provide more intuitive representation of obtained result.

Complementary analysis of the result from other methods.

Screening Preliminary identification of sensitive inputs. The methods are relatively easy to implement. The methods can be applied for

Non-quantitative analysis with many inputs.

Refined Complex analysis of sensitive inputs. Require an expertise and

relatively more resources Quantitative analysis.

Qualitative A heuristic characterization of parameters sensitivity. Require

relatively fewer model evaluations. Ranking input parameters. Quantitative Evaluation of the impact of model inputs to the variance in the

output. Require many model runs.

Deterministic or probabilistic analysis with few inputs.

(26)

26

Definition 3.1 The ANOVA-representation of function 𝑓(𝑥), with input space defined inside a unit n-dimensional cube [0,1]𝑛, is:

𝑓(𝑥) = 𝑓₀+ ∑ ∑ 𝑓_𝑖₁_…𝑖_𝑠(𝑥_𝑖₁, … , 𝑥_𝑖_𝑠)

𝑖1<⋯<𝑖𝑠 𝑛

𝑠=1

, (3.1)

where 𝑓₀ is a constant and the properties (3.2) and (3.3) are satisfied: 𝑓₀ = ∫ … 1 0 ∫ 𝑓(𝑥)𝑑𝑥 1 0 . (3.2) ∫ 𝑓_𝑖₁_…𝑖_𝑠(𝑥_𝑖₁, … , 𝑥_𝑖_𝑠) 1 0 𝑑𝑥_𝑖_𝑘 = 0, (3.3) where 1 ≤ 𝑘 ≤ 𝑠.

Since any input space can be transformed onto the unit hypercube, there is no loss of generality. Sobol proved that the decomposition (3.1) is unique (Sobol, 1990). The proof is reproduced in Appendix.

Example 3.1 (ANOVA-representation) Let us consider the function

𝑓(𝑥, 𝑦) = 𝑥𝑦, (3.4) where 𝑥 ∈ [0,1] and 𝑦 ∈ [0,1]. We will show how to find all terms of the ANOVA decomposition of the function 𝑓(𝑥, 𝑦).

𝑓₀ = ∫ ∫ 𝑓(𝑥, 𝑦)𝑑𝑥𝑑𝑦 1 0 1 0 = ∫ ∫ 𝑥𝑦𝑑𝑥𝑑𝑦 1 0 1 0 = (𝑥 2 2) 0 1 ∫ 𝑦𝑑𝑦 1 0 = (𝑥 2 2) 0 1 (𝑦 2 2) 0 1 =1 4. (3.5) 𝑔₁(𝑥) = ∫ 𝑥𝑦𝑑𝑦 = 1 0 𝑥 (𝑦 2 2) 0 1 =𝑥 2. (3.6) 𝑓₁(𝑥) = 𝑔₁(𝑥) − 𝑓₀ = 𝑥 2− 1 4. (3.7)

(27)

27 𝑔₂(𝑦) = ∫ 𝑥𝑦𝑑𝑥 = 1 0 (𝑥 2 2) 0 1 𝑦 =𝑦 2. (3.8) 𝑓₂(𝑦) = 𝑔₂(𝑦) − 𝑓₀ = 𝑦 2− 1 4. (3.9) 𝑓₁₂(𝑥, 𝑦) = 𝑓(𝑥, 𝑦) − 𝑓₀− 𝑓₁(𝑥) − 𝑓₂(𝑦) = 𝑥𝑦 −1 4− ( 𝑥 2− 1 4) − ( 𝑦 2− 1 4) = 𝑥𝑦 −𝑥 2− 𝑦 2+ 1 4. (3.10) The ANOVA-representation of the function 𝑓(𝑥, 𝑦) is

𝑓(𝑥, 𝑦) = 𝑓₀+ 𝑓₁(𝑥) + 𝑓₂(𝑦) + 𝑓₁₂(𝑥, 𝑦) = 1 4+ ( 𝑥 2− 1 4) + ( 𝑦 2− 1 4) + (𝑥𝑦 − 𝑥 2− 𝑦 2+ 1 4). (3.11) Indeed, by simplifying the equation (3.11) is the same as the initial function (3.4). The Sobol method applies the ANOVA representation of a function 𝑓(𝑥) in order to evaluate the Sobol sensitivity indices. To define the Sobol indices let our model described by the following function 𝑓: ℝ𝑛 → ℝ, 𝑿 ⟼ 𝑌 = 𝑓(𝑿), where 𝑓(∙) is the (unknown) model function integrable on [0,1]𝑛, 𝑿 = (𝑥₁, … , 𝑥_𝑛) are n independent random input parameters, and 𝑌 is the model response.

Since the terms in the ANOVA-representation are orthogonal, we can obtain the expected value of model response by

𝐸[𝑌] = ∫ … 1 0 ∫ 𝑓(𝑿)𝑑𝑿 1 0 = 𝑓₀. (3.12)

Moreover, for any point 𝑿 generated from n-dimension uniform distribution on [0,1]𝑛, the total variance of function 𝑓(𝑿) can be expressed by

𝑉 = 𝑉(𝑌) = ∫ … 1 0 ∫ 𝑓2(𝑿)𝑑𝑿 1 0 − 𝑓₀2. (3.13)

(28)

28

To define the partial variance, we use the property (3.3) to see that

(∫ … 1 0 ∫ 𝑓_𝑖₁_…𝑖_𝑠(𝑥_𝑖₁, … , 𝑥_𝑖_𝑠) 1 0 𝑑𝑥_𝑖₁_…𝑖_𝑠) 2 = 0. (3.14)

Definition 3.2 The partial variance corresponding to a term in (3.1) is given by

𝑉_𝑖₁_…𝑖_𝑠 = ∫ … 1 0 ∫ 𝑓_𝑖2₁_…𝑖_𝑠(𝑥_𝑖₁, … , 𝑥_𝑖_𝑠) 1 0 𝑑𝑥_𝑖₁… 𝑑𝑥_𝑖_𝑠. (3.15)

Then, squaring the decomposition (3.1), and taking into account the orthogonality of the terms, we obtain the equality

𝑉 = ∑ ∑ 𝑉_𝑖₁_…𝑖_𝑠. (3.16)

𝑖1<⋯<𝑖𝑠 𝑛

𝑠=1

The total deviation V determines the variation of the function 𝑓(𝑿), and the quantity 𝑉_𝑖₁_…𝑖_𝑠 characterize the impact of terms 𝑓_𝑖₁_…𝑖_𝑠 on the total variation.

Definition 3.3 Sobol global sensitivity indices are given by

𝑆_𝑖₁_…𝑖_𝑠 = 𝑉𝑖1…𝑖𝑠

𝑉 . (3.17)

It is clear that the value of an index 𝑆_𝑖₁_…𝑖_𝑠 is nonnegative, and less than one. In the thesis, the first order and total global sensitivity indices have been evaluated.

Definition 3.4 The first order Sobol sensitivity indices are the ratio

𝑆_𝑖 = 𝑉𝑖

𝑉, (3.18) with

(29)

29 𝑉_𝑖 = ∫ 𝑓_𝑖2(𝑥_𝑖)

1 0

𝑑𝑥_𝑖. (3.19)

In (Homma, et al., 1996), Homma et al introduced the total sensitivity indices.

Definition 3.5 The total sensitivity indices are defined as

𝑆_𝑖𝑡𝑜𝑡 = 𝑉𝑖

𝑡𝑜𝑡

𝑉 , (3.20) where 𝑉_𝑖𝑡𝑜𝑡 is the partial deviation, which includes the contribution of all the terms of the decomposition which depends on the variable 𝑥_𝑖 and can be obtained by

𝑉_𝑖𝑡𝑜𝑡 = 𝑉_𝑖 + ∑ 𝑉_𝑖,𝑗

𝑛 𝑗=1

𝑗≠𝑖

+ ⋯ + 𝑉_{1,…,𝑖,…,𝑛} (3.21)

Example 3.2 (Sobol indices) Let us find first-order and total sensitivity indices for the function (3.4). 𝑉 = ∫ ∫ (𝑥𝑦)2_{𝑑𝑥𝑑𝑦} 1 0 1 0 − 𝑓₀2 = (𝑥 3 3) 0 1 (𝑦 3 3) 0 1 − (1 4) 2 = 1 9− 1 16 = 7 144. (3.22) 𝑉₁ = ∫ 𝑓₁2(𝑥) 1 0 𝑑𝑥 = ∫ (𝑥 2− 1 4) 2 𝑑𝑥 = 1 0 (𝑥 3 12) 0 1 − (𝑥 2 8) 0 1 + ( 𝑥 16)0 1 = 1 48. (3.23) 𝑉₂ = ∫ 𝑓₂2(𝑦) 1 0 𝑑𝑦 = 1 48. (3.24) 𝑉₁₂ = ∫ ∫ 𝑓₁₂2(𝑥, 𝑦)𝑑𝑥𝑑𝑦 1 0 1 0 = ∫ ∫ (𝑥𝑦 −𝑥 2− 𝑦 2+ 1 4) 2 𝑑𝑥𝑑𝑦 1 0 1 0 = 1 144. (3.25) The first order Sobol sensitivity indices for 𝑥 and 𝑦 are equal to

𝑆₁ = 1 48 7 144 = 3 7. (3.26)

(30)

30 𝑆₂ = 3 7. (3.27) 𝑆₁₂ = 1 144 7 144 = 1 7. (3.28) 𝑆₁𝑡𝑜𝑡 = 𝑆₁+ 𝑆₁₂ =3 7+ 1 7 = 4 7. (3.29) 𝑆₂𝑡𝑜𝑡 = 𝑆₂+ 𝑆₁₂ = 4 7. (3.30) We can easily see that our computations are correct by sunning up the indices

𝑆₁+ 𝑆₂+ 𝑆₁₂ =3 7+ 3 7+ 1 7= 1. (3.31) The most informative extreme values of the sensitivity indices are:

A. 𝑆_𝑖 = 𝑆_𝑖𝑡𝑜𝑡 corresponds to the absence of the terms which depends on 𝑥_𝑖 along with other variables.

B. If 𝑆_𝑖 = 𝑆_𝑖𝑡𝑜𝑡 = 1, then the uncertainty in 𝑓(𝑿) depends only on the variable 𝑥_𝑖. C. If 𝑆_𝑖 = 𝑆_𝑖𝑡𝑜𝑡 = 0, then the uncertainty in 𝑓(𝑿) does not depend only on the variable

𝑥_𝑖.

D. If ∑𝑛_𝑖=1𝑆_𝑖 = 1, then 𝑓(𝑿) has additivity structure: 𝑓(𝑿) = 𝑓₀+ ∑𝑛_𝑖=1𝑓_𝑖(𝑥_𝑖).

Usually, Monte Carlo simulations are used to compute all required quantities (𝑓₀, 𝑉, 𝑉_𝑖₁_…𝑖_𝑠) in order to evaluate the sensitivity indices. In the thesis, a sample was produced using a quasi-Monte Carlo technique with the Sobol sequence. Details of this sampling approach are discussed in the next chapter.

(31)

31

4 Sampling methods

Traditionally, the methods presented in Chapter 2 and Chapter 3 use Monte Carlo (MC) simulations to estimate the corresponding values of interest. As depicted in Figure 4.1 (Becker, et al., 2015), MC methods are the appropriate techniques for models, which are not computationally expensive. In (Zijlema, 2010), it is noted that the simulations on unstructured meshes can increase the computational cost. In our case, the model SWAN requires approximately two minutes to run one simulation in stationary mode, and about 20 minutes in non-stationary mode for a 49 hours period on four processors. Therefore, we can assume that we deal with a computationally expensive model, and we must use approaches that are more advanced than crude MC methods.

Figure 4.1 – Methods for Sensitivity analysis (Becker, et al., 2015)

The MC approach is based on samples of pseudorandom numbers. In the thesis, the quasi-Monte Carlo method with Sobol sequences was used for uncertainty quantification and sensitivity analysis. First, we introduce the notions of random, pseudorandom and quasi-random numbers. Next, the basic concepts of the Sobol sequences are presented. Since any MC method is an approximation of the true solution,

(32)

32

error estimation of obtained result is an important issue, and we discuss this point in the last parts of this chapter.

One of the important steps in the design of the MC simulations is the generation of samples of input parameters for the multiple evaluations of the given function. Several most frequently used sample generation approaches are discussed in the following part.

4.1 Random and quasi-random numbers

A sequence of real numbers 𝛾₁, 𝛾₂… are called uniform random, if the variables are independent of each other, and uniformly distributed in the interval (0,1). Digital computers are used to produce pseudorandom sequences to simulate random numbers. In the thesis, the worlds “pseudorandom” and “random” are synonyms. The MC method with random sample has a convergence rate of the order 𝒪 ( 1

√𝑁). In order to obtain a faster

convergence, different variance reduction techniques exist, such as the Latin hypercube or orthogonal sampling. Hybrid-Monte Carlo methods also can be applied for high dimensional simulation (Ökten, 2001).

Quasi-random sequences are used instead of random numbers in order to obtain more uniform samples. In fact, quasi-random numbers are not random, but they are used for the same purpose and have some common properties with random sequences. In such a way, they have both the flexibility of pseudorandom numbers and the advantages of a regular grid. In addition, in some cases, they allows to obtain a better convergence rate in the MC integration due to its uniformity, and, consequently, higher accuracy. The MC methods with non-random sequences are called Quasi-Monte Carlo (QMC) methods.

In (Homma, et al., 1995), performances of the MC algorithm was examined with random numbers, the Latin hypercube sampling and the Sobol 𝐿𝑃_𝜏 sequences. Conclusions were made according to the mean square error of the obtained result, and it was shown that the 𝐿𝑃_𝜏 sequences have a better performance. Moreover, the crude MC and with Latin hypercube sampling perform approximately the same result.

(33)

33

Figure 4.2 shows graphically different sampling techniques, where we can see that the Latin hypercube and the Sobol sequences are more uniformly distributed than the random sample of points. The notion of uniformly distributed non-random points was introduced by Hermann Weyl in 1916.

Definition 4.1 A sequence of points in a unit hypercube 𝑥₁, 𝑥₂, … , 𝑥_𝑁 is uniformly distributed if for an arbitrary bounded Riemann-integrable function 𝑔(𝑥) (Sobol, 1998)

lim 𝑁→∞ 1 𝑁∑ 𝑔(𝑥𝑖) 𝑁 𝑖=1 = ∫ … 1 0 ∫ 𝑔(𝑥) 1 0 𝑑𝑥. (4.1)

Figure 4.2 – Two-dimensional scatter plots of the first 10 points of (a) random; (b) Latin hypercube; (c) Sobol sequences

4.2 Sobol

𝑳𝑷

𝝉

sequences

In 1967 Ilya Meerovich Sobol introduced a (t, s)-sequence (Sobol, 1967), which later was called the Sobol sequence. The sequences are generalized Niederreiter sequences with base equal to two (Niederreiter, 1992). The Sobol numbers were first example of digital sequences. Figure 4.3 illustrates the distribution of points in the Sobol sequences.

(34)

34

(35)

35

The (t, s)-sequence belongs to the class of the low-discrepancy sequences. In order to introduce the concept, let us start with the definition of term “discrepancy”.

Definition 4.2 Discrepancy is a quantitative measure for uniformity of a sample distribution. For the set of points 𝑃 = {𝑥₁, 𝑥₂, … , 𝑥_𝑁} with sample size equal to N, distributed in a n-dimensional unit hypercube, the discrepancy is

𝒟(𝐽; 𝑁) = 𝑠𝑢𝑝

𝐽

|𝐴(𝐽; 𝑁)

𝑁 − 𝑉𝑜𝑙(𝐽)|, (4.2) with the rectangular n-dimensional region

𝐽 = [a₁, 𝑏₁) × … × [a_n, 𝑏_𝑛) ⊆ [0,1]𝑛_{, (4.3)}

where 𝐴(𝐽; 𝑁) is the number of points in 𝐽 and 𝑉𝑜𝑙(𝐽) is the volume of 𝐽.

Example 4.1 (Discrepancy of a random sample) Consider the sample illustrated in Figure 4.4. The width of the small squares is equal to 0.4 and each of the square occupies 2/5 × 2/5 = 4/25 of the area of the unit square. In total, there are 25 randomly distributed points, hence, in the case of uniformly distributed sample, each of small square should contain 4 points. However, we see that one of the squares containing 7 points, and another one containing no points at all. Thus, the discrepancy of the set is at least 7/25 − 4/25 = 3/25.

(36)

36

In other words, discrepancy determines the maximum absolute difference between the fraction of the volume of a square and the fraction of the number of points inside the square. Since it is a difficult task to analytically compute the discrepancy of a sequence, the concept of star-discrepancy was introduced.

Definition 4.3 Star-discrepancy is the worst-case discrepancy, defined by 𝒟_𝑁∗ = max

𝐽 |𝒟(𝐽; 𝑁)|. (4.4)

The upper bound can be asymptotically estimated by

𝒟_𝑁∗ =𝐶(𝑛) log(𝑁)

𝑛

𝑁 , (4.5) where the constant 𝐶(𝑛) depends only on the specific sequence. For example, for the Sobol (t,s)-sequences the formula is

𝐶(𝑛) = 2

𝑡

𝑛! (log(2)𝑛₎. (4.6)

In our case, for the five-dimensional Sobol sequence the factor 𝐶(𝑛) is equal to 2.47 ∗ 10−2, according to (Niederreiter, 1992).

The upper bound value of the star-discrepancy (4.22) determines whether or not the sequence belongs to the class of low-discrepancy sequences.

The main property of low-discrepancy sequences is that as the length of the sequence becomes sufficiently large, the discrepancy has reduced to the rate, which is theoretically optimal.

In (Sobol, 1998), Sobol suggested criteria of uniformly distributed samples:

1. Better uniformity as the sample size increases.

2. Uniformity of the distribution for fairly small sample size. 3. Fast and simple algorithm for generation of points.

(37)

37

By following the procedure presented in (Bratley, et al., 1988), to generate the j-th component of points in a Sobol sequence, let us consider a primitive polynomial of degree d over the Galois field 𝔽₂, where 𝔽₂ = {0,1}

𝑃(𝑥) = 𝑥𝑑𝑗 _{+ 𝑎} 1 (𝑗) 𝑥𝑑𝑗−1 _{+ ⋯ + 𝑎} 𝑑_𝑗−1 (𝑗) 𝑥 + 1, (4.7) where the coefficients 𝑎₁(𝑗)… 𝑎_𝑑−1(𝑗) are either 0 or 1. These coefficients are used in order to define a sequence of positive numbers {𝑚_1,𝑗, 𝑚_2,𝑗… } by the following recurrence relation

𝑚_𝑖,𝑗 = 2𝑎₁(𝑗)𝑚_{𝑖−1,𝑗}⨁22𝑎₂(𝑗)𝑚_𝑖−2⨁ … ⨁ 2𝑑𝑗_𝑚

𝑖−𝑑𝑗,𝑗⨁𝑚𝑖−𝑑𝑗,𝑗, (4.8)

for 𝑑_𝑗 + 1 ≤ 𝑖, where ⨁ is the exclusive-or operator. The initial values 𝑚_1,𝑗, 𝑚_2,𝑗… , 𝑚_𝑑_𝑗_,𝑗 should be chosen such that each 𝑚_𝑖,𝑗, for 1 ≤ 𝑖 < 𝑑_𝑗 is a positive odd number less then 2𝑖. Now, let us define the direction numbers, {𝑐_𝑖} as

𝑐_𝑖,𝑗 = 𝑚𝑖,𝑗

2𝑖 , (4.9)

where 𝑐_𝑖_𝑗 is the binary representation of 𝑐_𝑖. The j-th component of i-th point is given by

𝑥_𝑖,𝑗 = 𝑏₁𝑐_1,𝑗⨁𝑏₂𝑐_2,𝑗⨁ …, (4.10) where (… 𝑏₂𝑏₁)₂ is the binary representation of i.

Example 4.2 (Generation of a Sobol sequence) Consider the primitive polynomial

𝑃(𝑥) = 𝑥3+ 𝑥 + 1. (4.11) The correspondent recurrence is

(38)

38

Let 𝑚_1,𝑗 = 1, 𝑚_2,𝑗 = 3, 𝑚_3,𝑗 = 7. By applying the recursion (4.8) we obtain 𝑚_4,𝑗 = 5, 𝑚_5,𝑗 = 7 and so on. This yields the direction numbers

𝑐_1,𝑗 = (0.1)₂, 𝑐_2,𝑗 = (0.11)₂,

𝑐_3,𝑗 = (0.111)₂, (4.13) 𝑐_4,𝑗 = (0.0101)₂,

𝑐_5,𝑗 = (0.00111)₂,

etc. Using (4.10) find the j-th component of the Sobol points

0 = (0)₂, 𝑥_0,𝑗 = 0, 1 = (1)₂, 𝑥_1,𝑗 = (0.1)₂ = 0.5, 2 = (10)₂, 𝑥_2,𝑗 = (0.11)₂ = 0.75, 3 = (11)₂, 𝑥_3,𝑗 = (0.1)₂⨁(0.11)₂ = (0.01)₂ = 0.25, (4.14) 4 = (100)₂, 𝑥_4,𝑗 = (0.111)₂ = 0.875, 5 = (100)₂, 𝑥_5,𝑗 = (0.1)₂⨁(0.111)₂ = (0.011)₂ = 0.375.

The first ten points of a ten dimensional sample obtained from the Sobol sequences are presented in Table 4.1.

(39)

39

4.3 Error estimation in Monte Carlo simulations

To estimate the error of the MC algorithm, we first review the Central Limit Theorem.

Theorem 4.4 (Central Limit Theorem) Let 𝑥₁, 𝑥₂, … , 𝑥_𝑁 be a sequence of independent identical distributed random variables of size 𝑁 with expected value 𝜇 = 𝐸(𝑥_𝑖) and finite variance 𝜎2 = 𝑉(𝑥_𝑖). Let 𝑋̅_𝑁 be the sample mean defined by

𝑋̅_𝑁 = 𝑥1+ 𝑥2 + ⋯ + 𝑥𝑁

𝑁 , (4.15) and 𝑍_𝑁 be the relation

𝑍_𝑁 = √𝑁𝑋̅𝑁 − 𝜇

𝜎 . (4.16) Therefore,

𝑍_𝑁 ⟹ 𝒩(0,1), (4.17) where ⟹ denotes the convergence in distribution.

In practice, the Central Limit Theorem helps to estimate so-called confidence interval. We note that the variable 𝑄_𝑁(𝑓) in equation (4.3) is an independent identically distributed random variable, and according to the Central Limit Theorem, the variable

√𝑁𝑄𝑁_𝜎−𝐼 has asymptotically a normal distribution.

Therefore, we can see that the approximation error 𝜀 = |𝑄_𝑁(𝑓) − 𝐼(𝑓)|, lies within the interval (−𝑧𝛼

2 𝜎 √𝑁, 𝑧𝛼2

𝜎

√𝑁) with probability 100 ∗ (1 − 𝛼)%, where the values of 𝛼 and

𝑧𝛼 2

can be taken from the standard normal table presented in most of the statistical

literature.

In (Sobol, 1990), Sobol suggested to use the probable error 𝛿_𝑁 = 0.6745 ∗ 𝜀_𝑁, where 𝜀_𝑁 = 𝜎

(40)

40

𝑃{|𝑄_𝑁(𝑓) − 𝐼(𝑓)| < 𝛿_𝑁|} ≈ 𝑃{|𝑄_𝑁(𝑓) − 𝐼(𝑓)| > 𝛿_𝑁|} ≈ 0.5. (4.18)

4.4 Error estimation in Quasi-Monte Carlo simulations

The exact upper bound of error of the QMC integration using the Sobol sequences can be approximated using the Koksma-Hlawka inequality:

|𝐼(𝑓) − 𝑄𝑁(𝑓)| ≤ 𝒟𝑁∗𝑉𝐻𝐾(𝑓), (4.19)

where 𝑉_𝐻𝐾(𝑓) is the variation of the function 𝑓(𝑥) in the sense of Hardy and Krause. For a one-dimension function 𝑓(𝑥), which first derivative is continuous, the variation has the form

𝑉_𝐻𝐾(𝑓) = ∫ |𝑑𝑓(𝑥) 𝑑𝑥 | 𝑑𝑥

1 0

. (4.20)

For higher dimensional functions, the variation may be expressed by an integral of partial derivatives, and, in general, it is not a trivial task to estimate the value of function variation. Nevertheless, from the inequality (4.19) together with (4.5) we can see that the upper bound of the approximation error in the result obtained by the QMC method

is 𝒪 (log(𝑁)𝑛

𝑁 ).

In (Archer, et al., 1997) and (Campolongo, et al., 2011), the bootstrap method for error estimation in the QMC simulations was used. We applied the same procedure for the results obtained in this research.

The bootstrap approach is based on resampling. Suppose we have received a sample {𝑥_𝑖𝑚} from QMC simulations. We resampled with replacement and repetitions the sample B times. For each of the new samples, the sensitivity indices 𝑆_𝑖 were estimated and we obtained the sample {𝑆_𝑖∗𝑏}_𝑏=1𝐵 , for 𝑖 = 1, … , 𝑛. Then, we applied the moment method to evaluate the 50% confidence intervals for 𝑆_𝑖 by

(41)

41 where 𝜀(𝑆̂_𝑖) = √ 1 𝐵 − 1∑(𝑆𝑖 ∗𝑏 _{− 𝑆} 𝑖∗⋆) 2 , 𝐵 𝑏=1 (4.22) with 𝑆_𝑖∗⋆ = 1 𝐵∑ 𝑆𝑖 ∗𝑏 𝐵 𝑏=1 . (4.23)

The result of uncertainty and sensitivity analysis obtained from the QMC simulations using Sobol low-discrepancy sequences is presented in the next chapter.

(42)

42

5 Results

The significant wave height was the model response variable under consideration in Uncertainty Quantification (UQ) and Sensitivity Analysis (SA). In order to estimate uncertainty and to illustrate its mutability with time, the model SWAN was run in non-stationary mode with half-hour time step. SA was implemented with respect to the results from runs in non-stationary and stationary modes. In UQ and SA, the result was obtained for the point corresponding to water gate ‘S-1’ of the Flood Prevention Complex depicted in Figure 1.1. In addition, the model outcome was obtained for each point of the unstructured mesh presented in Chapter 1, in order to obtain sensitivity maps.

5.1 Uncertainty Quantification

In the experiment for the UQ using Bayesian network approach, the following steps were executed:

1. Representing the model via a Bayesian network (Figure 5.1). 2. Assigning distribution to input parameters (Table 5.1).

3. Propagating uncertainty through the network using quasi-random samples. 4. Using the Quasi-Monte Carlo method, evaluating the integral

𝑓_𝑌

= ∫ 𝑓_𝑌(𝑦|𝑥₁𝑥₂𝑐₁𝑐₂𝑐₃)𝑓_𝑋₁(𝑥₁)𝑓_𝑋₂(𝑥₂)𝑓_𝐶₁(𝑐₁)𝑓_𝐶₂(𝑐₂)𝑓_𝐶₃(𝑐₃) 𝑑𝑥₁𝑑𝑥₂𝑑𝑐₁𝑑𝑐₂𝑑𝑐₃. (5.1) 5. Estimating the values of interest, such as mean, variance and standard deviation of the

model response.

The values of bottom level were represented on the unstructured mesh shown in Figure 1.4. In common practice, 10% variation is the suggested value of uncertainty in input parameters ( (Iooss, et al., 2015) and (Giap, et al., 2014)). Initially, uncertainty in the bottom level and wind parameters was 10%. Uncertainty in bathymetry was increased with growth of the depth value (see Table 5.1). In addition, we have executed an experiment with growing uncertainty in the wind velocity factor with time. The three

(43)

43

observed model coefficients had normal distribution. Since we did not calibrate these factors, we increased their uncertainty and set it equal to 30%.

Figure 5.1 – The model representation in the uncertainty study

Table 5.1 – Distribution of model inputs

Parameter Units Abbreviation Distribution Mean value Maximum deviation

Bottom level m BOTTOM Uniform -

10% for depth smaller than 50 meters;

40% for depth greater than 200 meters

Wind velocity m/s WIND Uniform -

10% for the experiment with constant uncertainty and the experiment in stationary mode;

58% at last period for the experiment with growing

uncertainty Bottom friction

coefficient - FRICTION Normal 0.0325 30% Whitecapping dissipation coefficient - cds2 Normal 0.000236 30% Wave steepness coefficient - stpm Normal 0.00302 30%

(44)

44

We executed two experiments with different settings for uncertainty in wind velocity: constant and increasing with time – each hour uncertainty was raised on 1% (Figure 5.2). Model uncertainty was estimated with respect to the result from 700 model runs for each of the experiments, with quasi-random samples, generated using the Sobol low-discrepancy sequence discussed in Chapter 4. As a result, we have seen that the standard deviation and uncertainty in the model response have increased significantly in the second experiment.

Figure 5.2 – Results of UQ

The measure of uncertainty (Figure 5.2, bottom plots) is equal to the ratio

𝑈 =𝜎

𝜇, (5.2) where 𝜎 is standard deviation, and 𝜇 is mean. The value of 𝑈 corresponds to the informational content of obtained forecast.

5.2 Sensitivity Analysis

The Sobol method for SA described in Chapter 3 was applied in order to determine the contribution of input parameters to model uncertainty. The goal was to answer the question: “Where is this uncertainty coming from?” (Saltelli, et al., 2010). We have

(45)

45

evaluated the importance of considered parameters in terms of the value of their Sobol indices: the more this value is close to 1.0, the more important this factor, and unimportant parameters have the value of their Sobol indices much less than 1.0.

SA was implemented according to the design suggested in (Saltelli, et al., 2010). A matrix of input parameters M of size (N, 2n) was generated, where N is the sample size, and n is the number of the input factors. The next step was to split the generated matrix

M into two matrices A and B, each of size (N, n). We defined the matrix 𝐴_𝐵(𝑖) equal to matrix A, except for the i-th column, which comes from matrix B, and the matrix 𝐵_𝐴(𝑖) was the same as matrix B, except for i-th column, obtained from matrix A. Model outputs for

all input values in 𝐴, 𝐵, 𝐴(𝑖)_𝐵 , and 𝐵_𝐴(𝑖) were denoted by 𝑓(𝐴), 𝑓(𝐵), 𝑓(𝐴_𝐵(𝑖)), and 𝑓(𝐵_𝐴(𝑖)), correspondingly.

The first order effects were estimated using Sobol’ indices given by

𝑆_𝑖=𝑉𝑋𝑖(𝐸𝑿~𝑖(𝑌|𝑋𝑖))

𝑉(𝑌) , (5.3) where 𝑋𝑖 denotes the i-th input parameter and 𝑿~𝑖is the matrix of all factors but 𝑋𝑖.

In order to compute the value of 𝑉𝑋_𝑖(𝐸_𝑿~𝑖(𝑌|𝑋𝑖)), we used the formula proposed in

(Sobol, et al., 2007): 𝑉_𝑋_𝑖(𝐸_𝑿~𝑖(𝑌|𝑋_𝑖)) ≃ 1 𝑁∑𝑓(𝐵)𝑗(𝑓(𝐴𝐵 (𝑖) )_𝑗− 𝑓(𝐴)_𝑗) 𝑁 𝑗=1 . (5.4)

The total sensitivity indices for i-th parameters were computed using the following formula:

𝑆_𝑇𝑖=𝐸𝑋~𝑖(𝑉𝑋𝑖(𝑌|𝑋~𝑖))

𝑉(𝑌) , (5.5) where the term 𝐸_𝑿~𝑖(𝑉𝑋_𝑖(𝑌|𝑋~𝑖)) was approximated by

𝐸_𝑋~𝑖(𝑉𝑋_𝑖(𝑌|𝑋~𝑖)) ≃ 1 𝑁∑𝑓(𝐴)𝑗(𝑓(𝐴)𝑗− 𝑓(𝐴𝐵 (𝑖) )_𝑗) 𝑁 𝑗=1 . (5.6)

(46)

46

As suggested in (Sobol, 1990), in order to decrease the approximation error, it is reasonable to pre-estimate the mean value of 𝑓(𝐴), and then subtract the value from each term in equations (5.4) and (5.6).

Sensitivity indices from the non-stationary mode experiment are plotted in Figure 5.3 with the reference to the 95% confidence intervals (upper plots).

(47)

47

As we can see, in the results from the experiment with constant uncertainty in wind parameter, there are much more oscillations in the values of sensitivity indices in comparison to the results from the experiment with growing uncertainty. According to Figure 5.4, we can explain such behaviour using the fact that importance of the wind velocity factor raised significantly in the second experiment, and this parameter dominated strongly at each point of time. Simultaneously, in the experiment with constant uncertainty, some parameters had effects that are more valuable at some points of time, depending on the value of significant wave height, and, consequently, their values fluctuated considerably.

Figure 5.4 – Scatter plot of sensitivity indices versus the value of significant wave height

The results of Sobol SA from the stationary mode experiment are shown in Figure 5.5 for wind velocity prediction for October 29, 2013 6 am (forecast for significant wave height was about 0.74 meters) and for 12 pm (significant wave height was predicted to be slightly above 1.5 meters). The total number of performed simulations was equal to 4200

(48)

48

for each time point. The error bars in Figure 5.5 are the bootstrap 50% confidence intervals, which were estimated using the approach, presented in Chapter 4.

In addition, the values of each parameter are given in Table 5.2. The value of 𝑆_𝑖/𝑆_𝑇𝑖 in Table 5.2 tells either uncertainty comes from the first order term of the decomposition (3.1), or from higher order terms. Since all of the parameters have a value of 𝑆_𝑖/𝑆_𝑇𝑖 grater or slightly below 0.9, we can conclude that uncertainty comes mostly from the first order terms of the factorization.

Table 5.2 –Results of Sobol SA Significant wave height (𝐻𝑠) Parameter 0.74 meters 1.5 meters 𝑆𝑖 𝛿(𝑆𝑖) 𝑆𝑇𝑖 𝛿(𝑆𝑇𝑖) 𝑆𝑖 𝑆𝑇𝑖 𝑆𝑖 𝛿(𝑆𝑖) 𝑆𝑇𝑖 𝛿(𝑆𝑇𝑖) 𝑆𝑖 𝑆𝑇𝑖 Bottom level 0.06 0.009 0.06 0.009 1.0 0.114 0.012 0.116 0.016 0.981 Wind velocity 0.613 0.032 0.62 0.022 0.989 0.588 0.031 0.594 0.023 0.99

Bottom friction coefficient 0.011 0.004 0.012 0.004 0.876 0.016 0.004 0.013 0.004 1.183

Whitecapping dissipation

coefficient 0.067 0.009 0.07 0.009 0.962 0.059 0.009 0.06 0.008 0.987

(49)

49

Figure 5.5 – First order and total sensitivity indices for five input parameters of the model SWAN

Figure 5.6 illustrates the convergence of first order and total sensitivity indices values to its average. The first 200 values of the sensitivity indices were treated as a warm up period and the average value was calculated using the last 400 values. These graphs show that with the sample size equal to 600, we have obtained the convergence to sample average.

A sensitivity map (Figure 5.7) contains the values of sensitivity indices calculated for each point of the computation grid. The sensitivity maps for bathymetry shows that uncertainty in the value of bottom level plays a significant role in some coastal areas. Because of this, we can conclude that the use of unstructured meshes for the representation of this factor can be treated as a very appropriate approach.

On the contrary, the significance of uncertainty in the value of wind velocity is approximately uniform, except some coastal regions. Therefore, this data can be presented on a regular mesh without loss of model response accuracy.

(50)

50

The sensitivity maps for the three observed coefficient similarly showed approximately uniform distribution of parameter effect with weak dependence on location.

Figure 5.6 – Convergence of first order and total sensitivity indices (SI) to their average value

(51)

51

Figure 5.7 – Sensitivity maps of first order Sobol indices for the region near the Flood Prevention Facility Complex

(52)

52

CONCLUSIONS

The results from the first experiment on Uncertainty Quantification show that uncertainty in the model response does not exceed 10% for the prediction of significant wave height larger than 0.5 meters. The results from the second experiment, where uncertainty in the wind parameter was growing with time, show the considerable increase of uncertainty in the model output and, as it was expected, the growth of importance of the wind speed parameter.

The sensitivity maps which we obtained for the bathymetry coefficient indicate the correlation between the location and the parameter importance. We can use flexibility of unstructured meshes and construct one such mesh so that its resolution is high in the areas where the output is more sensitive to uncertainty in this input, and decrease resolution for the rest of computational region.

Since the great importance of the wind velocity parameter was shown, in future work we will apply the study to the SWAN model coupling together with the HIRLAM model in order to find the precise value of uncertainty in this model input. The long-term goal is to build a real time prediction system, where uncertainty in the predicted value will be estimated and re-evaluated, whenever new wind velocity forecast is available.

(53)

53

APPENDIX

Theorem A.1. For any function 𝑓(𝑥) integrable on [0,1]𝑛, the decomposition of the form 𝑓(𝑥) = 𝑓₀+ ∑ ∑ 𝑓_𝑖₁_…𝑖_𝑠(𝑥_𝑖₁, … , 𝑥_𝑖_𝑠) 𝑖1<⋯<𝑖𝑠 , 𝑛 𝑠=1 (A. 1) with properties 𝑓₀ = ∫ … 1 0 ∫ 𝑓(𝑥)𝑑𝑥 1 0 . (A. 2) ∫ 𝑓_𝑖₁_…𝑖_𝑠(𝑥_𝑖₁, … , 𝑥_𝑖_𝑠) 1 0 𝑑𝑥_𝑖_𝑘 = 0, (A. 3) exists and it is unique.

Proof. Let us start by proving that the decomposition (A.1) of any integrable

function exists. First, consider the functions

𝑔_𝑖(𝑥_𝑖) = ∫ … 1 0 ∫ 𝑓(𝑥)𝑑𝑥_~𝑖 1 0 , (A. 4) where 𝑑𝑥_~𝑖 denotes the integration over all variables but xi . Let us define the first-order terms in (A.1) as following

𝑓_𝑖(𝑥_𝑖) = 𝑔_𝑖(𝑥_𝑖) − 𝑓₀. (A. 5) In order to show that such definition of the first-order terms satisfies the property (A.3), we integrate (A.5) over xi.

∫ 𝑓_𝑖(𝑥𝑖)𝑑𝑥𝑖 = 1 0 ∫ (𝑔_𝑖(𝑥𝑖) − 𝑓0)𝑑𝑥𝑖 = 1 0 ∫ 𝑔_𝑖(𝑥𝑖)𝑑𝑥𝑖 − 𝑓0 = 1 0 𝑓₀− 𝑓₀ = 0. (A. 6)

In this way, we have shown that (A.3) is satisfied, and that we have built the first-order terms of the decomposition (A.1).