Efficient uncertainty quantification of
a fully nonlinear and dispersive
water wave model with random inputs
Abstract
A major challenge in nextgeneration industrial applications is to improve numerical analysis by quantifying uncertainties in predictions. In this work we present a formulation of a fully nonlinear and dispersive potential flow water wave model with random inputs for the probabilistic description of the evolution of waves. The model is analyzed using random sampling techniques and nonintrusive methods based on generalized Polynomial Chaos (PC). These methods allow to accurately and efficiently estimate the probability distribution of the solution and require only the computation of the solution in different points in the parameter space, allowing for the reuse of existing simulation software. The choice of the applied methods is driven by the number of uncertain input parameters and by the fact that finding the solution of the considered model is computationally intensive. We revisit experimental benchmarks often used for validation of deterministic water wave models. Based on numerical experiments and assumed uncertainties in boundary data, our analysis reveals that some of the known discrepancies from deterministic simulation in comparison with experimental measurements could be partially explained by the variability in the model input. We finally present a synthetic experiment studying the variance based sensitivity of the wave load on an offshore structure to a number of input uncertainties. In the numerical examples presented the PC methods have exhibited fast convergence, suggesting that the problem is amenable to being analyzed with such methods.
Keywords:
uncertainty quantification generalized polynomial chaos sensitivity analysis highperformance computing free surface water waves.Msc:
65C99 76B151 Introduction
In coastal and offshore engineering it is important to design maritime structures that can withstand critical failures due to waveinduced loadings. The most extreme wave inducedloadings can be estimated from direct measurements, laboratory experiments and simulationbased tools, which can account for the wave kinematics sufficiently accurately. It is still common to predict wave kinematics using numerical tools that have been validated by a single or few deterministic simulations and compared to idealized physical experiments, e.g., in wave tanks. These validation procedures, and at a greater scale real field simulations, are subject to a number of uncertainties that could lead to unexpected results. The study of the influence of these uncertainties on the resulting engineering analysis and decisions requires ultimately a shift from a deterministic approach to a probabilistic approach Wojtkiewicz2001. Such engineering analysis are useful for risk management aimed at reducing risk in design and operations.
The research field of Uncertainty Quantification (UQ) includes all the techniques used to rigorously study uncertainties entering a systems and their influence on its dynamics. These techniques deliver tolerance intervals and probability distributions of system outputs, denoted Quantities of Interest (QoIs). Upon the knowledge of a deterministic model describing the dynamical system, denoted the forward model, uncertainty quantification can be split in four steps:

identification of sources of uncertainty and QoIs,

quantification of uncertainty sources by means of probability distributions,

uncertainty propagation through the system,

sensitivity analysis.
Uncertainties in coastal and offshore engineering are either related to weather conditions or to structural/bathymetry characteristics. In the first case, the weather conditions are commonly grouped into a number of sea states characterized by a number of parameters, which determine particular probability distributions associated to waves, wind, currents, sea level and ice characteristics. Transitions between different sea states are modeled by stationary processes governing these parameters. The characterization of these uncertainties – step (ab) – requires extensive measurements and the current state of the art is presented in BitnerGregersen2014; BitnerGregersen2014a.
A significant part of existing works on the propagation of uncertainty and sensitivity analysis of water waves use the Shallow Water Equations (SWE) as forward model, thus addressing mostly tidal waves, where vertical velocities are negligible. The focus is usually on the characterization of extreme responses or, in other words, the probability of failure of certain safety conditions. In Naess2012 such analysis are performed using the MC method. In Ge2008 the Monte Carlo (MC) and the Polynomial Chaos (PC) method in the Galerkin and collocation form are compared when applied on the SWE modeling the propagation of a wave over a submerged hump. In Ge2011 the PC method in its collocation form is used for the study of uncertainties in floodhazard mapping. In Liu2009 random sampling methods, sparse grid, PC in Galerkin and collocation form, and a novel quadrature technique called Compound Uncorrelated Dimension (CUD) quadrature were compared when applied on the SWE modeling flood prediction under an uncertain river bed topography and characteristics. In Ricchiuto2014 a combination of nonintrusive (collocation) PC and ANOVA decomposition was used for the propagation and sensitivity analysis of the uncertain parameters entering the SWE modeling the runup of waves. Unlike the preceding works, Yildirim2015 studies the influence of uncertainties on the phaseaveraged equation, which is suitable for slowly varying wave fields, e.g. ocean waves in deep water. Yildirim2015 considers uncertainties entering the source term, the boundary conditions and the current field, adopting PC and ANOVA decomposition approaches.
1.1 Paper contributions
We consider a fully nonlinear and dispersive potential flow model, which allows the study of the phase resolved propagation of water waves over varying bathymetry. The model has traditionally been consider computationally very costly, however, recent progresses in the design of scalable linear solvers has made the model much more tractable for improved analysis of predictions via uncertainty quantification methods. UQ analysis require an overall computational cost which often scales badly, due to the curse of dimensionality, with the computational cost of the forward model. Thus, efficient numerical solvers for the forward model and efficient UQ techniques must be used. We employ state of the art software EngsigKarupEtAl2008; EngsigKarupEtAl2011; GlimbergEtAl2011; EGNL13; LEDN13 to obtain efficient and accurate simulations of the forward model, while we turn to recently developed techniques based on polynomial chaos Xiu09fastnumerical; Maitre2010 to perform uncertainty quantification.
These techniques will be applied to classical benchmarks, such as BB94; DutykhClamond2013, where different sources of uncertainty and QoIs will be investigated. Due to the lack of data, some assumptions will be made about the probability distributions of the sources of uncertainty – step (b), that in hydrodynamics simulations are commonly inlet/outlet conditions (boundary conditions), bathymetry data and structural positions (geometry). In experimental settings, all these uncertainties can be classified as epistemic Kiureghian2009, because they can in principle be reduced either by better measurements and/or by more precise experimental designs.
We will propagate the uncertainties through the dynamical system – step (c) – using traditional sampling techniques, such as Monte Carlo methods, and modern techniques based on generalized Polynomial Chaos (gPC) Xiu09fastnumerical; Maitre2010. Nonintrusive approaches such as stochastic collocation and sparse grid will be preferred to intrusive approaches, due to the ability of the former of reusing existing code, avoiding the need for reengineering existing software.
For one of the analyzed test cases we will identify the sources of uncertainty to which the QoI is most sensitive – step (d). Techniques addressing this problem exist and we will focus on the variance based method of Sobol’ Sobol1993; Rabitz2000; Chan2000; Sudret2008; Crestaux2009; Alexanderian2012, where the sensitivity of a QoI to a certain input uncertainty is expressed by the amount of variance of the QoI which is due to such input uncertainty. The identification of the least influential inputs allows the refinement of the model, where the uncertainties on these inputs are disregarded. The method will be used to quantify the sensitivity to a number of input uncertainties of the wave load on an offshore structure. The sensitivities can be used to refine the forward model, disregarding uncertainties on parameters that do not influence significantly the QoIs.
1.2 Paper organization
The paper is organized as follows. In section 2 we introduce the governing equation for the deterministic description of dispersive and nonlinear water waves based on potential theory. In section 3 we describe how a model with random inputs can be formulated and parametrized. Section 4 presents the theory of random sampling methods and of polynomial chaos methods for the forward propagation of uncertainty, and the description of the method of Sobol’ for sensitivity analysis. In section LABEL:sec:UQonWaterWaves, the effect of parametric uncertainty in bathymetry and wave input are studied and numerical experiments are compared for different approaches. Section LABEL:sec:Conclusions contains concluding remarks.
2 Mathematical formulation of the forward model
We consider unsteady water waves described by a potential model for threedimensional fully nonlinear and dispersive free surface flows under the influence of gravity. The flow is assumed inviscid and irrotational. It can, without simplifications, be used for short and long wave propagation in both shallow and deep water where viscous and rotational effects are negligible. The sea bed is assumed variable and impermeable.
We introduce a Cartesian coordinate system with the horizontal and the vertical dimensions, where the coordinate points upwards. The functions and describe respectively the depth of the sea bed and the free surface. The still water level is given by .
2.1 The deterministic potential flow model
The evolution of water waves over an arbitrary sea bed are described by the kinematic and dynamic free surface boundary conditions^{1}^{1}1The gravitational acceleration constant, , is set to be 9.81 .
(1a)  
(1b) 
where . We will consider waves in a spatial domain (fluid volume), and a time domain with final time . For the fluid volume, a Laplace problem defines the scalar velocity potential
(2a)  
(2b)  
(2c) 
Using a classical transformation
(3) 
the Laplace problem can be written as
(4a)  
(4b)  
(4c) 
where
(5a)  
(5b)  
(5c) 
The relation between the scalar velocity potential function and velocity field is
(6) 
The governing equations can be solved in the setting of a numerical wave tank and is then subject to initial and boundary conditions
(7) 
where wave generation and absorption is done using a line relaxation method LD83. A complete derivation of the equations are given in EGNL13. These model equations can be solved numerically using flexibleorder finite differences EngsigKarupEtAl2008; EngsigKarupEtAl2011 and the massively parallel implementation LEDN13 enables fast hydrodynamics computations EGNL13. A fast solver is a prerequisite for enabling UQ within acceptable time frames for realistic engineering applications.
2.2 Calculation of loads on structures
The estimation of the pressure in the water column is needed for load calculations. We start from the momentum equation
(8) 
where is the local pressure and is the fluid density. By integration of the pressure in the vertical direction and assuming that the flow is inviscid and irrotational it is possible to derive the pressure equation
(9) 
The remaining integral term can be estimated numerically, e.g. by using a sufficiently accurate numerical quadrature rule.
For the estimation of structural forces, we can integrate the expression for the pressure given in Eq. (9) in the vertical direction, to obtain the force as
(10) 
which for a vertical structure takes the form
(11a)  
(11b)  
(11c) 
The integral terms can be estimated numerically.
3 The model with random inputs
Despite the effort towards the emulation of idealistic phenomena, experimental settings suffer from unavoidable epistemic uncertainties. In the context of water waves, these uncertainties characterize the experimental geometry as well as the boundary conditions. The parameters describing these uncertainties are collected into a dimensional vector of random variables defined on a probability space describing such uncertainties, where is the space of elementary events, is a field and is a probability measure^{2}^{2}2For a treatment of measure theoretic probability theory see e.g. Billingsley2008. In the following we will avoid the direct use of the probability space by the introduction of the image space , where is the Borel algebra over and is the cumulative distribution function (CDF) of . This will allow the evaluation of integrals of integrable functions of by noting that .
The reformulation with random inputs of (1) and (2) leads to a system of PDEs for the unknowns and , which are now random fields:
(12a)  
(12b) 
3.1 Uncertain bathymetries, random fields and KarhunenLoève Expansion
The bathymetry function describing stillwater depth in the considered domain can be uncertain. In laboratory experiments, this can be due to manufacturing errors, whereas in the real settings this is often due to the lack of precise knowledge of the bottom topography or the unknown action of sedimentation. Since this uncertainty is spatially varying, it must be treated as a random field. We will treat this random field using the KarhunenLoève (KL) expansion 1977probability; Schwab2006 which is a useful parametrization technique when the field is characterized by a certain amount of spatial correlation.
Let be a spatially varying random field over a spatial domain with mean and covariance function . Then the bathymetry function can be parametrized as an infinite series
(13) 
where the convergence is in , are random variables, , and are the solutions of the generalized eigenvalue problem
(14) 
If is a Gaussian random field, then .
For practical computations (13) is truncated at a desired order . It is easy to check how much of the variance of the original random field is retained by such approximation, using that
where the orthogonality of is exploited. There are several options regarding the correlation kernel . All these are problem dependent and an appropriate characterization of the random field has to be performed prior to the construction of the KLexpansion. In this work, we will use the exponential covariance kernel
(15) 
where is the correlation length. By this choice is then an OrnsteinUhlenbeck Uhlenbeck1930 random field. The eigenvalue problem (14) is solved numerically using a spectral discretization of the integral operator. Figure 1 shows realizations of the KLexpansions of a 1D random field for the submerged bar experiment considered in section LABEL:sec:submergerbar with exponential covariance kernel and zero mean for different correlation lengths . The total variance represented by the KLexpansions is fixed to (the total variance of with exponential covariance kernel is ). In figure (a)a and (b)b, fields with different correlation lengths are illustrated. Shorter correlation lengths determine a slower decay of the expansion coefficients in (13) and thus an expansion retaining an higher number of terms is required to express higher local variability.
4 Uncertainty Quantification
We are interested in studying the propagation of uncertainties through the dynamical system (12). To reduce the notation used, let . We describe the results in terms of the probability distribution and/or the first moments, e.g., mean and variance, of :
(16)  
(17) 
where is the Probability Density Function (PDF) of the random vector .
In this work we will focus exclusively on nonintrusive methods, which require a minimal development effort. In particular the existing solvers are considered as black boxes and the nonintrusive methods need only to be wrapped around them. On the contrary, intrusive methods require the development of new solvers based on mixed discretization of the stochastic and the spatial models. These methods are usually better in dynamically adapting to timedependent problems Cheng2013; Cheng2013a; Boyaval2010; Venturi2006; Sapsis2009 but their implementation is often very demanding – sometimes prohibitive – for existing customized solvers.
4.1 Pseudorandom sampling methods
Among the existent UQ techniques, pseudorandom sampling methods are the most widely used. The most notable of these techniques, is the Monte Carlo (MC) method. It is based on the law of large numbers which states that given the random vector and the functional ,
(18) 
for . In the definition above is a ensemble of samples drawn from the probability distribution of and a.s. stands for almost surely implying convergence in probability. The probabilistic error of these approximations is reduced asymptotically as . In spite of this slow convergence rate, MC method is widely used due to its robustness, ease of use and to the fact that it does not suffer the curse of dimensionality, i.e., its convergence rate is independent from . Improvements of the MC method have been proposed in order to cover more uniformly the parameter space, obtaining improved convergence rates.
The Latin Hypercube Sampling (LHS) method Mckay2000 divides the parameter space in equiprobable bins along each dimension and samples are taken such that each bin contains only one sample. This provides a better convergence rate in the average cases, even if the worst case convergence rate remains .
The Quasi Monte Carlo (QMC) method Morokoff1995 generates low discrepancy sequences of points such that the domain is uniformly covered. The convergence rate of QMC is improved to , but the dimension of the parameter space becomes relevant.
4.2 Deterministic sampling methods
In the following we will handle functions with finite variance, i.e. belonging to the weighted space defined as
(19) 
with inner product and norm defined as
(20) 
For many standard onedimensional distributions with density , we can find the set of univariate polynomials that form an orthonormal basis for Xiu:2002:WPC:587159.587325; Xiu2010. For instance, orthogonal basis for the space of random variables with Beta, Gaussian or Gamma distributions are the Jacobi, the probabilists’ Hermite and the Laguerre polynomials respectively. If the distribution is not standard, but admits a probability density function, then one can still use GramSchmidt orthogonalization to create suitable polynomials, or use the Stieltjes procedure to compute the recursion coefficients of the polynomials orthogonal with respect to such measure (see Gautschi2004; Gautschi1994).
For random vectors of mutually independent random variables with densities and basis , it holds that . Thus, a possible set of basis functions for is given by the set of multivariate polynomials where is a multiindex and
(21) 
This construction includes basis functions and is more accurate than the required order . An alternative set of basis is given by , where . For this set of basis N+dN