Bivariate Kernel Density Estimation for Meta-ocean Contour Lines of Extreme Sea States

-,-

(a.State Key Laboratory of Ocean Engineering;b.School of Naval Architecture,Ocean and Civil Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

Abstract:A novel environmental contour lines(meta-ocean contour lines)method was proposed based on measured ocean wave data at two offshore sites. The use of bivariate kernel density estimation with biased cross-validation bandwidth selection was proposed for implementing the novel environmental contour lines method.Comparison of the environmental contours obtained by using the proposed novel method with those obtained by using the Clayton copula transformation method clearly substantiated the effectiveness and superiority of the proposed novel method. The research results demonstrate that the proposed novel method can be utilized as an effective tool for predicting the long-term extreme dynamic responses of ocean engineering structures.

Key words:environmental contours;bivariate kernel density estimation;biased cross-validation bandwidth selection;Clayton copula transformation

0 Introduction

In order to successfully design an offshore structure (an offshore rig, an offshore wind turbine or a wave energy converter)，it is of vital importance to predict the severity of sea states that the offshore structure may encounter during its operational life (e.g.50 years).Current practice for designing an offshore structure often involves applying nonlinear time domain numerical simulations to predict the response of the structure to extreme sea states. Offshore structure designers often use the following risk assessment procedures:(1)Represent the sea area of interest using buoy observations or hindcast simulations of an acceptable period,typically about 20 years;(2)Use extreme value theory and models to extrapolate to more extreme values of waves from a shorter period of record;(3) Generate meta-ocean contour lines (environmental contour lines) involving significant wave heightHSand spectral peak periodTP(or energy periodTe)to determine extreme sea states;(4)Identify one or more sea states to describe an individual or group of waves to use as an input to a numerical or physical model simulation; (5) Perform simulations or testing to predict the offshore structure response in these conditions.

The determination of the extreme sea states is the most important step in the afore-mentioned design procedures. These meta-ocean contour lines predict future extreme sea states based on a finite amount of existing data gathered from buoys measurements or hindcast simulations. Metaocean contour lines must predict the envelope of sea states likely to occur in some return period,which is often longer than the period of record. A good meta-ocean contour should neither overpredict or under-predict the bounds of this envelope, as this could result in costly over-design or failure.

In order to successfully design and build an offshore structure, it is imperative to be able to forecast the long-term sea conditions that the offshore structure will have to endure in a given location. Creating meta-ocean contours based on hindcast buoy data was first proposed by Haver[1]in 1987 and is still a frequently used method of accomplishing this.Meta-ocean contour lines take ordered pairs of observed sea state data and determine an area of observations in which exceeding them could potentially elicit extreme structural responses in a specified return period (for higher combinations ofHSandTe).The study in this paper uses two commonly-used variables—significant wave height (HS) and energy period (Te). According to the theory in the field of ocean engineering,a wave’s height(HS)will be related to/dependent upon a variable relating to its width (Te).In order to generate a meta-ocean contour line, a joint distribution of the two variables must be established.The dependency of the two variables makes modeling a joint distribution far more challenging.Since the original meta-ocean contour line method was first proposed, many additional contour methods have been developed[2-3,5-8]. Where these methods all primarily differ is the way in which they attempt to capture the dependency between the two variables.

If the joint probability distribution ofHSandTeis available,the Inverse-First-Order Reliability Method(IFORM)with the Rosenblatt transformation[9]can be applied to derive a meta-ocean contour line of extreme sea states[10]. Because the Rosenblatt transformation is applied, the joint probability distribution ofHSandTeis usually described by means of a set of conditional distributions.In the process of deriving the meta-ocean contour line, Rosenblatt transformation relates the environmental variables(HSandTe)in the physical parameter space with independent standard normal variables (u1,u2) in the so-called standard normal space. This Rosenblatt transformation-based metaocean contour line method was used by Haver and Winterstein[5]to assess the accidental tether loads in a tension leg platform and used by Berg[11]on position mooring. However, oftentimes, the 50-year meta-ocean contour lines created by using this method with the Rosenblatt transformation fail to cover the measured data taken over a relatively short time period (8～20 years).The shortcomings of the contours created by using this method with the Rosenblatt transformation obviously demonstrate that this method needs further modification and enhancement for generating more realistic meta-ocean contours.

In the field of ocean engineering,a marginal probability distribution is usually modeled for significant wave heightHS,and then the probability distribution of the spectral peak periodTP(or energy periodTe)is modeled conditional on the significant wave heightHS.As the problem dimension increases, the amount of data required for estimating the conditional distribution parameters will increase substantially. Therefore, it may become very difficult to apply these models when the problem dimension increases. Consequently, in the real world ocean engineering problems it is more common to estimate the marginal distributions of the random variables and their correlation structures based on the statistical data.

Manuel et al[12]developed meta-ocean contour lines by considering environmental variables with a Clayton copula. They derived analytical expressions for the meta-ocean contour lines in terms of the marginal distributions of the environment variables and their rank-related correlation coefficients. Their proposed method can be conveniently implemented and provides a straightforward way to easily produce meta-ocean contours for more than two environment variables. However, their method has not been verified by extensive measured ocean wave data in different offshore sites.

With the motivation to overcome the shortcomings of the aforementioned methods, in this study, we will propose the use of bivariate kernel density estimation with rigorous bandwidth selection for generating meta-ocean contour lines.The meta-ocean contour lines at two offshore sites obtained by using our proposed new method will be compared with those obtained by using the Clayton copula transformation method, and the superiority and effectiveness of our proposed new method will be fully and clearly substantiated.

This paper begins in Chapter 1 by elucidating the theories behind the traditional methods and our proposed new method for generating meta-ocean contour lines.It continues in Chapter 2 by applying our proposed new method for generating meta-ocean contour lines based on measured ocean wave data at two offshore sites. These results are compared with those created by using the Clayton copula transformation method,with concluding remarks summarized in Chapter 3.

1 Meta-ocean contour line(environmental contour)methods

1.1 Clayton copula-based environmental contours

The meta-ocean contour line method can be applied for an offshore site if the joint probability distribution function for the significant wave heightHSand the energy periodTeis available in the form of a joint modelFHSTe( )h,t. However, the modelling of the joint probability distribution is typically very difficult, and therefore, in the field of ocean engineering, a marginal probability distribution is usually modeled for the significant wave heightHS, and then the probability distribution of the energy periodTeis modeled conditional on the significant wave heightHS. In the meantime, the correlation structures between the random variables should also be estimated. When Kendall’sτ,a rank-related correlation coefficient,has been estimated and is available along with marginal probability distributions for the significant wave heightHSand the energy periodTe, a Clayton-copula transformation may be employed to define meta-ocean contour lines.

In statistics,the Kendall rank correlation coefficient,commonly referred to as Kendall′s tau coefficient (after the Greek letterτ), is a statistic used to measure the ordinal association between two measured quantities. It is a measure of rank correlation: the similarity of the orderings of the data when ranked by each of the quantities. Intuitively, the Kendall correlation between two variables will be high when observations have a similar (or identical for a correlation of 1) rank (i.e. relative position label of the observations within the variable: 1st, 2nd, 3rd, etc.) between the two variables,and low when observations have a dissimilar (or fully different for a correlation of -1)rank between the two variables.

In probability theory and statistics, a copula is a multivariate cumulative distribution function for which the marginal probability distribution of each variable is uniform. Copulas are used to describe the dependence between random variables.Sklar’s theorem states that any multivariate joint distribution can be written in terms of univariate marginal distribution functions and a copula which describes the dependence structure between the variables.

In this paper, without loss of generality, a specific copula family (i.e., Clayton copula) is chosen to derive copula-based meta-ocean contour lines, which are then compared with those based on our proposed novel method, The calculation results from these two methods are consistent with the same met-ocean data for the selected two offshore sites.

Contour lines corresponding to a constant annual exceedance probability can be obtained by transforming the model to a space consisting of independent, standard normal variables (u1,u2) by using the copula transformation. In the standard normal space, the two-dimensional meta-ocean contour line associated with a return periodTr(in years)is a circle of radiusβthat can be calculated as follows:

whereΦ( )is the standard normal distribution function andTsis the sea state duration expressed in hours.

The copula transformation requires estimating the Kendall’s τ coefficients between the physical variablesHSandTeaccording to Eq.(1).

The copula transformation first maps the random variables (HS,Te) onto the standard normal variables(u1,u2)as follows:

Then from Eq.(9)we have:

The Clayton copula distribution function is expressed as follows:

The relationship between the parameterθof the Clayton copula and the Kendall’sτcoefficient is as follows:

The Kendall’sτcoefficients can be calculated based on the measured values of the physical variables(HS,Te),and the values of the parametersθcan then be calculated by utilizing Eq.(15).

Finally, in the Clayton copula transformation, the primary random variable, the significant wave heightHS,is mapped tou1directly:

whereγis a randomly-generated angle with a uniform distribution in[0,2π].

In the Clayton copula transformation, the secondary random variable, the energy periodTe, is mapped to the normal random variableu2as follows:

By transforming the circles in the normal space back to the physical parameter space using the above two equations,we have finally obtained the meta-ocean contour lines.

1.2 Bivariate kernel density estimation with rigorous bandwidth selection

The use of the methods such as those listed in Section 2.1 requires a predefined probability model for the considered environmental variables. These methods are frequently constrained by their underlying structural assumptions. Oftentimes the observations of sea state variables have great variations, which requires a more flexible approach not constrained by the predetermined assumptions of the contour behavior.In this study,we propose the use of bivariate kernel density estimation with rigorous bandwidth selection as a novel method for generating theHS-Temeta-ocean contour lines.

In statistics,kernel density estimation is a non-parametric way to estimate the probability density function of a random variable.Kernel density estimation is a fundamental data smoothing problem where inferences about the population are made, based on a finite data sample. The general equation for a multivariate kernel density estimation ford-dimensional datax1, …,xdgiven multivariate data setX1,…,Xnis as follows[13]:

whereKis the kernel—a non-negative function,andh＞0 is a smoothing parameter called the bandwidth parameter. Intuitively one wants to choosehas small as the data will allow; however, there is always a trade-off between the bias of the estimator and its variance.The bandwidth of the kernel is a free parameter which exhibits a strong influence on the resulting estimate.The most common optimality criterion used to select the bandwidth parameterhis the mean integrated squared error(MISE):

wherefis the generally-unknown real density function. Under weak assumptions onfandKwe have the following relationship[14]:

2 Calculation examples and discussions

2.1 The calculation example based on the measured data at National Data Buoy Center Station 42039

Our first calculation example is based on the measured data at National Data Buoy Center(NDBC)Station 42039(see Fig.1).The NDBC 42039 buoy is located at a sea site 115 nautical miles from Pensacola,Florida,U.S.A..(https://www.ndbc.noaa.gov/station_page.php?station=42039).The water depth at this measuring site is 270 m . Fig.2 shows the position (largest yellow mark) of the National Data Buoy Center Station 42039. In this study, the methodologies described in Chapter 1 will be applied to 182 607 hourly observations ofHSandTetaken from January 1, 1996 to December 31,2018 at NDBC 42039.

Fig.1 National Data Buoy Center Station 42039

Fig.2 Position(yellow mark)of the National Data Buoy Center Station 42039

Tab.1 shows a part of the derived significant wave heights and energy periods at NDBC 42039.

Tab.1 Derived significant wave heights and energy periods at NDBC 42039

Fig.3 shows the 50-year extreme sea state contours created by using the Clayton copula transformation method and the new approach presented in this paper for NDBC42039. In Fig.3 the red‘dots’represent the measured data.In Fig.3 the black dashed curve represents the 50-year extreme sea state contour created by using the Clayton copula transformation method (i.e.using Eqs.(1)-(20)implemented in the WDRT toolbox. The WDRT (WEC Design Response Toolbox) software was developed by Sandia National Laboratories and the National Renewable Energy Laboratory (NREL) to provide extreme response analysis tools, specifically for design analysis of ocean structures such as wave energy converters (WECs). It is expected that, given a period of record on the order of about 23 years,the extreme contour for a return period of 50 years should include all of the observed data.However,the black dashed contour created in Fig.3 using the Clayton copula transformation method obviously does not follow the shape of trends present within the measured data. This black contour is not successful at fitting smaller values ofHSandTe,but also does not allow for coverage in the area of the input space in which theHSandTevalues are relatively high.

Fig.3 50-year extreme sea state contours created by the Clayton copula method and the proposed new method

In Fig.3 the green dashed curve represents the 50-year extreme sea state contour created by using the bivariate kernel density estimation method implemented in the WDRT toolbox. However,the WDRT toolbox does not have the capability of calculating the bandwidth parameterhthat is required (in Eq.(21)) during the kernel density estimation process. Therefore, in our study, we have imported the measured NDBC42039 data into MATLAB in the form of a Hierarchical Data Format(i.e. NDBC42039.h5). Then we have calculated a bandwidth parameterh(=0.149 7) for one variate(HS) in MATLAB by using the biased cross-validation bandwidth selection method based on implementing Eqs.(21)～(27) in Chapter 1. Similarly, we have also calculated another bandwidth parameterh(=0.352 8)for another variate(Te)in MATLAB by using the biased cross-validation bandwidth selection method based on implementing Eqs.(21)～(27) in Chapter 1. After obtaining the values of these two bandwidth parameters, we subsequently exported them into the WDRT toolbox. By running the WDRT code we finally obtained the green dashed curve in Fig.3 representing the 50-year extreme sea state contour.

This green dashed contour created using the bivariate kernel density estimation with rigorous bandwidth selection is very successful at fitting both smaller values ofHSandTeand higher values ofHSandTe. As is expected, given a period of record on the order of about 23 years, the green dashed extreme contour for a return period of 50 years has included all of the observed data (except one data point with a very largeTevalue).Here it should be noted that this missed data point with a very largeTevalue has a relatively small significant wave height.It is therefore practically insignificant because the long-term extreme structural dynamic responses are generally caused by a sea state with a very high significant wave height. Obviously,the green dashed contour in Fig.3 follows the shape of trends present within the measured data.Most noticeably,this green dashed contour allows for coverage in the area of the input space in which theHSvalues are very high.

2.2 Calculation example based on the measured data at NDBC Station 42035

Our second calculation example is based on the measured data at NDBC Station 42035 (see Fig.4). The NDBC 42035 buoy is located at an offshore site 22 nautical miles from Galveston, Texas, U.S.A. (https://www.ndbc.noaa.gov/station_page.php?station=42035). The water depth at this measuring site is 16.2 m. Fig.5 shows the position (the largest yellow mark) of the NDBC Station 42035. In this study, the methodologies described in Chapter 1 will be applied to 188 748 hourly observations ofHSandTetaken from January 1,1996 to December 31,2018 at NDBC 42035.

Fig.4 National Data Buoy Center Station 42035

Fig.5 Position(the largest yellow mark)of the NDBC Station 42035

Tab.2 shows a part of the derived significant wave heights and energy periods at NDBC 42035.

Tab.2 Derived significant wave heights and energy periods at NDBC 42035

Fig.6 shows the 50-year extreme sea state contours created by the Clayton copula transformation method and the new method presented in this paper for NDBC42035. In Fig.6, the red‘dots’represent the measured data, and the black dashed curve represents the 50-year extreme sea state contour created by using the Clayton copula transformation method (i.e. using Eqs.(1)～(20)) implemented in the WDRT toolbox. It is expected that, given a period of record on the order of about 23 years, the extreme contour for a return period of 50 years should include all of the observed data.However, the black dashed contour in Fig.6 created using the Clayton copula transformation method obviously does not follow the shape of trends present within the measured data. This black contour is not successful at fitting smaller values ofHSandTe. This black contour does not allow for coverage in the area of the input space in which theHSvalues are relatively high. This black contour also does not allow for coverage in the area of the input space in which theTevalues are relatively high.

Fig.6 50-year extreme sea state contours created by the Clayton copula method and the proposed new method

In Fig.6, the green dashed curve represents the 50-year extreme sea state contour created by using the bivariate kernel density estimation method implemented in the WDRT toolbox. However,the WDRT toolbox does not have the capability of calculating the bandwidth parameterhthat is required(in Eq.(21))during the kernel density estimation process.Therefore,in our study,we imported the measured NDBC42035 data into MATLAB in the form of a hierarchical data format (i.e.NDBC42035.h5). Then we calculated a bandwidth parameterh(=0.145 5) for one variate (HS) in MATLAB by using the biased cross-validation bandwidth selection method based on implementing Eqs.(21)～(27)in Chapter 1.Similarly,we also calculated another bandwidth parameterh(=0.291 3)for another variate (Te)in MATLAB by using the biased cross-validation bandwidth selection method based on implementing Eqs.(21)～(27)in Chapter 1.After obtaining the values of these two bandwidth parameters, we subsequently exported them into the WDRT toolbox. By running the WDRT code, we have finally obtained the green dashed curve in Fig.6 representing the 50-year extreme sea state contour.

This green dashed contour created using the bivariate kernel density estimation with rigorous bandwidth selection is very successful at fitting both smaller values ofHSandTeand higher values ofHSandTe. As is expected, given a period of record on the order of about 23 years, the green dashed extreme contour for a return period of 50 years has included all of the observed data (except two data points with very largeTevalues).Here it should be noted that these missed data points with very largeTevalues have very small significant wave height values. They are therefore practically very insignificant because the long-term extreme structural dynamic responses are generally caused by a sea state with a very high significant wave height. Obviously, the green dashed contour in Fig.6 follows the shape of trends present within the measured data. Most noticeably, this green dashed contour allows for coverage in the area of the input space in which theHSvalues are high.

3 Concluding remarks

With the motivation to overcome the shortcomings of the traditional methods, in this study, we proposed the use of bivariate kernel density estimation with biased cross-validation bandwidth selection for generating meta-ocean contour lines. The meta-ocean contours at two offshore sites obtained by using the proposed novel method were compared with those obtained by using the Clayton copula transformation method, and the effectiveness and superiority of the proposed novel method were clearly substantiated. The research results in this paper demonstrate that the proposed novel method can be utilized as an effective tool for predicting the long-term extreme dynamic responses of ocean engineering structures.