Deriving environmental contours from highest density regions

Environmental contours are an established method in probabilistic engineering design, especially in ocean engineering. The contours help engineers to select the environmental states which are appropriate for structural design calculations. Defining an environmental contour means enclosing a region in the variable space which corresponds to a certain return period. However, there are multiple definitions and methods to calculate an environmental contour for a given return period. Here, we analyze the established approaches and present a new concept which we call highest density contour (HDC). We define this environmental contour to enclose the highest density region (HDR) of a given probability density. This region occupies the smallest possible volume in the variable space among all regions with the same included probability, which is advantageous for engineering design. We perform the calculations using a numerical grid to discretize the original variable space into a finite number of grid cells. Each cell's probability is estimated and used for numerical integration. The proposed method can be applied to any number of dimensions, i.e. number of different variables in the joint probability model. To put the highest density contour method in context, we compare it to the established inverse first-order reliability method (IFORM) and show that for common probability distributions the two methods yield similarly shaped contours. In multimodal probability distributions, however, where IFORM leads to contours which are dificult to interpret, the presented method still generates clearly defined contours.

U Random variable in standard normal space [-] u Realization of the random variable in standard normal space [-] X Random variable in original space [-] x Realization of the random variable in original space [-] z

Purpose of environmental contours
Engineers have to design any marine structure in such a way that it is able to withstand the loads induced by the environment.As the environment, i.e. wind, waves and currents, continually change and cannot be predicted for long periods of time, the environment is often modeled stochastically by defining probability density functions, f (x j ).Then, the structure is designed to withstand all but some extremely rare environmental states, e.g.all waves with significant wave heights, H s , less than a threshold, h s , with a cumulative probability or exceedance probability of α, i.e.P r(H s ≤ h s ) = 1 − α or P r(H s > h s ) = α.In general notation for any random variable, X 1 , there exists a threshold, x 1 , which fulfills The exceedance probability, α, corresponds to a certain recurrence or return period, T , which describes the average time period between two consecutive environmental states above the threshold, x 1 .The threshold is called return value.For example, to comply with standards a marine structure such as an offshore wind turbine is required to withstand significant wave heights, H s , with a return period, T , of 50 years [20].Often, however, structural safety depends not only on one variable, but on the occurrence of combinations of p variables, {X j } p j=1 .When two variables are of importance, e.g.significant wave height, H s , and spectral peak period, T p , a joint probability density function can be defined and an environmental contour can be calculated which encloses the subset (or region) of environmental states that the structure has to be designed for.Here, we call this region design region (Fig. 1).Often the most critical structural response is associated with very high or low values of environmental variables, i.e. with environmental conditions located at the boundary of the design region.Consequently, standards allow engineers to calculate structural responses for a limited set of environmental design conditions along the contour instead of requiring engineering calculations based on a high number of possible variable combinations spread over the complete design region [8].If there are more than two variables the concept of environmental contours leads to environmental surfaces (3 variables) or environmental manifolds (> 3 variables).Here, for simplicity we also refer to these as environmental contours.

Different definitions and methods
As there are different mathematical definitions for environmental contours one has to further specify which kind of environmental contour is being calculated.Different concepts of environmental contours lead to different design loads and consequently to different structural responses [1].Originally, environmental contours arose from the concept of return values in univariate probability density functions which are calculated based on one-sided exceedance over threshold (Fig. 2a).Consequently, a logical definition for an environmental contour is (i) constant one-sided exceedance in all directions of the p-dimensional variable space, P r(X 1 > x 1 , X 2 > x 2 , ..., X p > x p ) = α.The bottom panel in Fig. 2a shows the contour for the twodimensional joint distribution of X 1 and X 2 .However, for design purposes not only the highest values of a variable can be of interest, but also the lowest.For example, when designing an offshore structure, low values of the peak period, T p , have to be considered as the structure's natural frequencies can be either higher or lower than the average peak period.Consequently, another possible definition for an environmental contour is (ii) two-sided exceedance over threshold (Fig. 2b; e.g.[21]).A third possibility is to define an environmental contour to have (iii) constant probability density, f m , along its path enclosing the most likely environmental states (Fig. 2c).In this case a T -year return period means that on average every T years an environmental state with a probability density less than f m occurs.In the broader statistics literature the variable region enclosed by such a contour is called a highest density region (HDR) [19].Although HDRs are a logical concept for environmental contours, yet no author has strictly followed this definition.The design curve introduced by Haver [14] is a related concept since it is a line of constant probability density, but only one-sided exceedance is considered.The constant probability density approach described by Det Norske Veritas [8] does define a fully closed contour of constant probability density.However, it is designed in such a way that it is unclear how much probability is enclosed by the contour.Instead the contour's probability density, f , is chosen to be the joint probability density of the (x 1 , x 2 )variable combination with x 1 = return value based on the marginal x 1 -distribution and x 2 = an associated x 2 value (Fig. 3c).Leira [23], however, has indeed used the HDR definition but only after a transformation of the original variables into standard normal space.When transforming the contour back to the original variable space the constant probability density is not preserved.Here we will compute contours strictly following the HDR definition.Besides these different definitions of types of environmental contours there exist different methods to calculate a given type of environmental contour.The traditional and probably most used approach is the so-called inverse first-order reliability method (IFORM) [34,15].It is a standard design practice for a wide range of marine engineering applications where extreme sea states are of interest [8].These are for example ships [11], offshore wind turbines [20], floating structures [9] or wave energy converters [7,10].Using IFORM one defines the contour in standard normal space, U j , instead of the original environmental variable space, X j .Thus, one first defines a circle with a radius, β, in the U -space (Fig. 3a).The radius corresponds to the return period and increases with longer periods.Then one transforms the points along the circle to the original variable space leading to the environmental contour.This transformation is done via the inverse Rosenblatt transformation [27].As its name implies IFORM is a reliability method and is based on the idea that the exceedance region approximates the failure region, F, of a structure (and the exceedance probability, α, approximates the structure's failure probability, P f ; see [25]).Contours based on IFORM are widely used and have been published e.g. by Saranyasoontorn and Manuel [28], Leira [23], Baarholm et al. [2], Li et al. [24], Myers et al. [26], Valamanesh et al. [30], Eckert-Gallup et al. [13].
Huseby et al. [17], however, pointed out that the Rosenblatt transformation introduces errors as failure probabilities, P f , can be underestimated or overestimated on a case by case basis.Therefore, they introduced an alternative    (a) Inverse first-order reliability method (IFORM) [34].
The contour is defined as a circle in standard normal variables, U j .The points along the circle have to be transformed to the original environmental variables, X j .(b) Huseby et al. [17] Monte Carlo contour.The contour is computed with the environmental variables, X j , directly.(c) Constant probability density contour described by Det Norske Veritas [8].By its definition it is unclear how much probability is enclosed by the contour.(d) Jonathan et al. [21] constant exceedance probability contour.The calculation can be done in a general set of variables, M j , e.g. in the X-or U -space.In comparison to the other methods a different definition for the exceedance region is used (compare shaded areas).
method to calculate environmental contours in the original variable space.Following their method, one first carries out a Monte Carlo simulation to generate a high number of sea states based on a given joint probability distribution model.Then one chooses an angle, θ, defining a line (in two dimensions, p = 2) and varies its position such that it divides the variable space into one halfspace containing the majority of data points and the other halfspace containing the data points representing the exceedance probability, α × n (with n being the total number of simulated environmental states, Fig. 3b).By iterating this procedure over a finite number of angles, θ ∈ [0, 360), the resulting lines can be connected to an environmental contour.This new approach has been picked up in several recent publications, e.g. to compare the approach to the traditional IFORM method [33], to compare different statistical models [31] or to decrease the required process time [18].While the Monte Carlo method overcomes the problems caused by the Rosenblatt transformation it requires the simulation of environmental states which is computationally more expensive than the simple IFORM calculations.Further, by its definition the method cannot generate concave contours.Jonathan et al. [21] define and calculate environmental contours yet differently.Using clear mathematical notation they find a contour with constant exceedance probability, P r(X 1 > x 1 , X 2 > x 2 ) = α (notation for two dimensions, p = 2).Thus, instead of finding halfspaces which are tangential to the contour, their exceedance regions have finite boundaries for each variable leading to outwards radiating rectangles in a two-dimensional Cartesian coordinate system (Fig. 3d ).Consequently, in contrast to IFORM and the Monte Carlo approach the method does not try to match the exceedance region with the failure region and thus separates the concept of an environmental contour from a structure's failure function.Following this method one first chooses a reference point, r M 0 .Then one defines a line which passes through that point at an angle, θ, to the abscissa.Lastly one finds the position along the line which satisfies P r = α.Repeating this procedure over a full circle, θ ∈ [0, 360), one finds the environmental contour.The method can be applied in any variable space, M j , e.g. in original variables, X j , or standard normal variables, U j .Further, besides fully closed contours, one-sided exceedance is also considered by the authors.One can combine the method with using modern conditional extreme models [16] as demonstrated by Jonathan et al. [22,21].The method disconnects the environmental state statistics from any particular structural problem which makes it a more general approach to define a T -year set of environmental states for any further use of these data.However, like the reliability methods, it defines multiple α-exceedance regions in the variable space of a single probability model.While in reliability methods the idea is that one of these multiple exceedance regions overlaps with the failure region this is not the case with the P r(X 1 > x 1 , X 2 > x 2 ) = α definition.Thus, if a contour is defined independently of the concept of failure regions, it seems more meaningful to define α to be the probability of a single region (instead of having multiple regions with α probability content).
Motivated by the individual advantages the described contour calculation methods have, here, we introduce contours enclosing highest probability density regions which we compute using numerical integration.We continue the idea introduced by Jonathan et al. [21] of decoupling the exceedance region from the structure's failure region, but go one step further and do not define any kind of outwards radiating exceedance region.Instead, we choose to find a contour which encloses the most likely environmental states which together make up a defined probability of 1 − α.The proposed method allows us to define the contour in the original variable space and can be used for any number p of dimensions.By discretizing the variable space into a finite number of grid cells and using numerical integration techniques any probability distribution can be evaluated, e.g common parametric sea state joint probability distributions [32], nonparametric models [12] or extreme value models which can have discontinuities at the threshold [29].Similar as being done e.g. in computational fluid dynamics [5] we demonstrate that with a sufficiently small grid cell size the solution is grid independent.

Data
In order to compare our environmental contour approach to similar methods we use the 3-hour sea state model presented by Vanem and Bitner-Gregersen [32].They use a fitted joint model for significant wave height, H s , and zero-upcrossing period, T z .Based on their model environmental contours have been calculated using both the traditional IFORM method [32] and the newer Monte Carlo method [17].The joint model was derived from one particular location in the ERA-Interim data set [6].Significant wave height, H s , is modeled as a 3-parameter Weibull distribution with the parameters α W (scale), β W (shape) and γ W (location): Based on a least squares fit the parameters are α W = 2.776, β W = 1.471 and γ W = 0.8888 [32].
The zero-upcrossing period, T z , is modeled to follow a log-normal distribution, LN : The distribution's parameters, µ Hs and σ Hs , are conditional on the significant wave height, H s , and are modeled as 3-parameter functions: In this case they are estimated to be a 1 = 0.1000, a 2 = 1.489, Multiplying the marginal distribution of the significant wave height, H s , and the conditional distribution of the zero-upcrossing period, T z , we can calculate the joint distribution: Since the data represent 3-hour sea states, exeedance probability, α, for a T -year return period is calculated as 3. Highest density contour (HDC)

Analytical definition
Our goal is to find a contour, C, of constant probability density, f m , which encloses a probability of 1 − α, i.e.: This contour, C, encloses the highest density region, R. Therefore we call C highest density contour (HDC).A highest density region fulfills two main properties: (i) the probability density of every point inside is at least as large as the probability density of any point outside and (ii) for a given probability content the region occupies the smallest possible volume in the variable space [4].There is no general analytic solution to find the HDR or HDC, i.e. solving for C or R in Eq. 8.
HDRs, however, can be computed based on numerical integration approaches [35] or Monte Carlo techniques [19].Environmental contours involve very low α values and are usually based on low-dimensional probability models.Thus, we choose numerical integration over Monte Carlo simulation to compute the highest density contour, C.However, if a probability model, which incorporates many environmental variables (high p value), is evaluated numerical integration might become infeasible and Monte Carlo approaches should be used.Here, we use numerical integration and start by discretizing the probability density space into a finite number of equally sized grid cells.In the next section we will evaluate the two-dimensional case, but in the appendix the equations for p dimensions are given.

Numerical integration approach in two dimensions
The two-dimensional probability space is discretized in K × L grid cells with a constant size of ∆x 1 × ∆x 2 (Fig. 4).Each grid cell's center point, (x 1k , x 2l ), is used as the reference position of the cell.Then, based on the cumulative distribution function, F X1 , the cell-averaged probability density in the first dimension, fX1 , is calculated using central difference: The cell-averaged probability density in the second dimension, fX2|X1 , is calculated similarly: While fX1 is the true cell-averaged probability density in the first dimension, in the second dimension, fX2 is a c c e p t e d v e r s i o n approximated since the dependence of F X2|X1 upon x 1 within the grid cell is not accounted for.Instead we fix x 1 to the value at the grid cell center, x 1 = x 1l , and therefore assume F X2|X1 to be constant from x 1l − 0.5∆x 1 to x 1l + 0.5∆x 1 .
Multiplying the two individual probability densities yields the cell-averaged joint probability density, f : Now we can compute the probability that an event with a minimum probability density of f m occurs, i.e. we calculate the probability content enclosed by a HDC of f m probability density.This probability, F (f m ), is calculated by summing up the probabilities of all cells which have a probability density greater than or equal f m (Fig. 4): If the joint probability density function is unimodal the grid cells which fulfill f ≥ f m make up a single contiguous area.The boundary of this area is a contour which encloses a probability of F .Using the function F (f m ) we can consequently find a contour with a given exeedance probability, α, of interest by finding the corresponding minimum probability density, f m : Solving this equation is a root finding problem of a monotonically decreasing function ( F (f m ) − 1 + α = 0).We solve the equation using Matlab's (version R2015b, The MathWorks, USA) fzero function which iteratively finds the root of a nonlinear function.All grid cells fulfilling f ≥ f m then approximate the HDR, R(f m ), and the grid cells making up the boundary of the HDR approximate the HDC, C(f m ).

Properties of the highest density contours
As done in previous work based on the described joint probability model [17,32] we compute the 1-, 10-and 25year environmental contours (Fig. 5).The corresponding exceedance probabilities are α 1 = 3.42 × 10 −4 , α 10 = 3.42 × 10 −5 and α 25 = 1.37 × 10 −5 respectively.The computed HDCs have constant probability densities of f m1 = 4.4 × 10 −5 (1-year), f m10 = 4.3 × 10 −6 (10-year) and f m25 = 1.7 × 10 −6 (25-year).Fig. 6a shows how the enclosed probability, F , monotonically decreases with increasing f m until it reaches F = 0. Since the probability functions we use here (Weibull and log-normal) are unbounded, F asymptotically approaches 1 as f m approaches 0. Fig. 6b presents the maximum H s -and T z -values along a contour of constant f m -probability density (H s , T z ).Longer return periods, T , lead to smaller f m -values and consequently to bigger contours with higher H s and T z values.
As discretization in general is sensitive to step size we evaluate the contour's robustness with respect to grid cell size ∆x 1 = ∆H s , ∆x 2 = ∆T z .We analyze how minimum probability density, f m , changes with grid cell size.In all three tested return periods (1-, 10-and 25-year contour) minimum probability density, f m , is roughly constant at small cell sizes and starts to fluctuate with increasing cell size indicating a grid-independent solution can be reached (Fig. 7a).Oversized grid cells can lead to minimum probability density being half or double than the converged minimum probability density (Fig. 7b).For the given probability model we find that convergence is reached at a grid cell size of H s = 0.05 m and T z = 0.05 s.There, deviation to the smallest tested grid cell size is less than 1 %, 0.99 < f * m < 1.01, with f * m being minimum probability density, f m , normalized by the converged f m value (Fig. 7c).

Comparison with IFORM and Monte Carlo contours
For comparison we further compute environmental contours using IFORM based on the same probability model.The highest density contours have similar shapes as the contours calculated with IFORM and the Monte Carlo method (Fig. 8c,d ).However, we define a HDC to enclose a probability of 1 − α while an IFORM contour and a Huseby et al. [17] Monte Carlo contour each enclose a probability less than 1 − α since by their definitions multiple regions outside the contour have a probability of α (Fig. 8a).Consequently, the HDC's dimensions in terms of H s and T z are bigger in comparison.However, for a fairer comparison we can inflate an IFORM contour and find  the T -year contour which encloses exactly 1 − α probability.Leira [23] showed that this can be done by utilizing the inverted Rayleigh distribution (for two dimensions).The author calls these contours equi-shape contours.Here, we find that such a 25-year equi-shape contour corresponds to a 308.8-yearIFORM contour.The contour's shape and size is roughly similar to the 25-year HDC.These similarities suggest that the 308.8-yearIFORM contour has approximately constant probability density, f m25 , along the contour.
To visualize a typical data set, we Monte Carlo simulate 25 years of 3-hour sea states (n = 73050; gray dots in Fig. 8c).In this particular data set one data point exceeds the HDC while there are multiple data points exceeding the 25-year IFORM contour.The different contour dimensions can also be expressed in terms of maximum H s -and T z -values along the contour (H s , T z ).While Huseby et al. [17] report 25-year maxima of H s 25 = 14.66 m and T z 25 = 13.68 s for the Monte Carlo contour, here we find H s 25 = 16.79 m and T z 25 = 14.64 s for the HDC and H s 25 = 15.23 m and T z 25 = 13.96s for the IFORM contour (Fig. 8d ).Thus, the HDC H s 25 value is 10.2 % higher than the IFORM value and 14.5 % higher than the Monte Carlo method value.Consequently, from an engineering design point of view the HDC is the most conservative method of the three considered.
This does not only apply to the considered distribution, but is a generic property based on the different definitions of these three contours.The IFORM and Monte Carlo contours are defined to contain the return value of the marginal distribution as their highest variable value, i.e.H s 25 = H s,25 (Fig. 8a).On the other hand, a HDC is defined to enclose 1 − α probability.Since it does not contain all H s -T z sea states fulfilling H s < H s,25 (which together would make up 1 − α probability) it must contain some sea states with H s > H s, 25 .
By the HDC's definition of an enclosed probability of 1 − α, in a random 25-year data set the probability that at least one data point exceeds a 25-year contour is about 63.2 %, 1−(1−α 25 ) n ≈ 0.632 with n = 25×365.25×24/3= 73050.Here, exceedance precisely means that this sea state realization is anywhere outside the region enclosed by the contour, R(f m ).Such a sea state occurs on average every 25 years.This simple and clear interpretation is why we have chosen the definition of constant probability density and a probability of 1−α, i.e. defining the contour to enclose the highest density region.We believe that this definition offers an intuitive and meaningful concept for a T -year environmental contour in the engineering design process.If an engineer designs a structure to withstand all sea states inside a T -year contour, the structure will be designed for the most likely (extreme) sea states which are expected to occur in T years.Then on average every T years a sea state will occur which the structure is not designed for.
Alternative concepts with multiple α-exceedance regions (see Fig. 3a,b) are based on the idea of known failure regions in the context of structural reliability methods (see [25]).IFORM assumes that a structure's failure surface (or limit state surface) has a convex shape.It defines the α-halfspace exceedance regions in its particular way because in that case the true failure surface can be linearized such that the variable space is separated by a straight line at an angle θ into a survival region and a failure region (in two dimensions).Then, this failure region overlaps with IFORM's exceedance region.It has the failure probability P f = α and the survival region the survival probability 1 − P f .Here, however, we completely separate the idea of describing the environmental conditions from any particular structural problem.Thus, we do not intend to align the α-probability exceedance region with a particular failure region.
As described IFORM leads to a contour which encloses less than 1 − α probability and consequently results in less conservative design conditions compared to a HDC.If the structural design, which is developed based on these environmental conditions, has a convex failure surface, the theoretical precondition of IFORM is met.Then in comparison, a HDC can be seen as overly conservative.Thus, if the designer knows that a structure responds with a convex failure surface choosing an IFORM contour is advantageous in the sense that it yields less conservative but still safe design conditions.
While many structures respond with a convex failure surface this precondition for IFORM connects the environmental contour to a certain class of structures.The shape of the failure surface might be unknown beforehand and 1-10      Sketches showing expected differences in contour size due to different definitions.Some contours are defined in such a way that the maximum value along the contour, H s 25 , is equal the return value of the marginal distribution, H s25 , (middle).The highest density contour (HDC), however, is defined to enclose 1 − α probability and thus has a maximum value along the contour which is higher than the return value of the marginal distribution, H s 25 > H s25 , (right).(b) Sketch illustrating an IFORM contour and possible failure regions of a linear system of three components, F 1 , F 2 , F 3 .Since the contour contains less than 1 − α probability the system's failure probability can be greater than α.(c) A total of n = 73050 sea state data points have been Monte Carlo simulated representing a 25-year data set (scatter plot).The 25-year HDC (solid line) and the 25-year IFORM contour (short dashes) have similar shapes, but as expected the HDC is bigger.The 308.8-year IFORM contour or 25-year equi-shape contour (long dashes; [23]) encloses the same amount of probability as the 25-year HDC.(d) Comparison of maximum values along the contour, H s and T z .As expected by the different definitions, the HDC has the highest maximum significant wave height, H s , and maximum zero-upcrossing period, T z .a c c e p t e d v e r s i o n only becomes apparent during the design process.If it turns out that the failure surface is non-convex and therefore violates IFORM's precondition the designer would need to go one step backwards and define new design conditions by inflating the IFORM contour.By not making use of the properties of possible structural responses the HDC is more conservative, but also more general in its application.It would avoid the need of the described iteration loop in the design process.
Further, a highest density contour is advantageous in the design process of a structural problem of a system consisting of multiple components.Consider a series structure consisting of z different components with z different failure functions.In a series structure a failure of one component results in failure of the system [3].Suppose that each component fulfills IFORM's precondition of having a convex failure surface.Nevertheless, the probability contained by the union of all z failure regions, F 1 ∪ F 2 ... ∪ F z , could exceed α (Fig. 8b).In that case it would be expected that frequenter than every T years an environmental state occurs which leads to failure of some of the components and consequently failure of the system.If an environmental contour containing 1 − α probability were used to design the components, on the other hand, by definition the system's probability of failure would be less than α.Consequently, the system would be expected to survive longer than T years.
A similar example could be given for a single component with multiple failure modes.The three failure regions shown in Fig. 8b would then correspond to different failure modes and the same conclusions as for the series structure could be drawn.These two examples explain why IFORM is primarily aimed at assessing the reliability of one component failing in one particular failure mode.A highest density contour, on the other hand, could be used in these two cases without worrying that any assumptions might be violated.

Bimodal mixture model
Highest density contours can be computed based on any probability distribution.The used definition of constant probability density along the contour, f m , can lead to multiple enclosed subregions for a given return period, T , if the probability distribution is multimodal.Here, we demonstrate this by extending the joint H s -T z -probability distribution by a mixture model for the zero-upcrossing period, T z .We use the H s -T z environmental variables although we are aware that such a H s -T z distribution might be physically unrealistic.However, for simplicity we build upon the previously described sea state model instead of setting up a new case with a different set of environmental variables.Thus, we keep the log-normal distribution term, LN ( µ Hs , σ 2 Hs ), from Eqs. 3-5 and mix it with a normal distribution, N (µ 2 , σ 2  2 ): Similar to the parameters µ Hs and σ Hs we define the mixture coefficient, p Hs , to be conditional on significant wave height, H s .Using an exponential decay function, we let the normal distribution term, N (µ 2 , σ 2 2 ), fade out at high significant wave height, H s : We design two mixture models.For the first model we create a normal distribution, N , such that its probability density blends smoothly into the log-normal distribution, LN , by using a mean value of µ 2 = 10 s and standard deviation of σ 2 = 2 s (model 1 ).For the second model we design a normal distribution which has much less density overlap by using a mean value of µ 2 = 15 s and standard deviation of σ 2 = 0.5 s (model 2 ).For both models we compute the 25-year HDC as well as the 25-year IFORM contour.In model 1 the HDC and IFORM contour have similar shapes.Both have a concave path at high T z -values and as expected the HDC is bigger in size (Fig. 9a).In contrast, model 2 has two distinct probability density maxima which lead to different shapes for the IFORM and HDC.While the HDC encloses two separated subregions the IFORM contour encloses a single contiguous region (Fig. 9b).This single region contains sea states with much lower probability densities than the conservative HDC as by its definition IFORM can only enclose one single contiguous region.Consequently, in this example an engineer who designs a structure to withstand all loads caused by the environmental states along this 25-year contour would design the structure to withstand some environmental states which are expected to occur extremely rarely.Therefore, possible structural designs which are limited by these environmental states would not be considered which could lead to bad design, either from a cost or engineering perspective.
The apparent difference in shape between the two contours is interesting since it visually demonstrates that the IFORM contour does not have constant probability along its path and consequently does not enclose the most likely environmental states.Strictly, this should not be expected anyway, but since it is roughly true for many ordinary sea state models, one might intuitively interpret an IFORM contour that way.By IFORM's definition the contour has two properties in the U -space: (i) constant probability density along its path and (ii) α-probability halfspaces separated by lines which are tangent to the contour (Fig. 3a).Interestingly, for many sea state probability models these two properties roughly translate to the X-space.Here, we demonstrate the rough persistence of the constant probability density property for the unmodified sea state model since in this case IFORM and HDCs have similar shapes (Fig. 8b).Rough persistence of the α-halfspace property, on the other hand, has been shown by Huseby et al. [17] who computed Monte Carlo contours which are defined by enforcing the α-halfspace property in the original variable space (Fig. 3b).These Monte Carlo contours have been reported to have similar shapes as the IFORM contours.Thus, based on experience an engineer might intuitively interpret a typical IFORM contour to have roughly constant probability density and α-halfspace exceedance probability in the original variable space.This interpretation would not hold true for the multimodal model 2, however.In addition to clearly not having constant probability density it also does not roughly have α-halfspace exceedance in the original variable space since the contour is concave.Not having any meaningful properties in the original variable space, raises the question how to intuitively interpret an IFORM contour in such a case.In contrast, the presented highest density contour with its constant probability density, f m , along the contour and its enclosure of a probability of 1−α offers a clear interpretation for any probability distribution.

Summary and conclusions
In this work we present environmental contours which enclose regions of highest probability density.A highest density contour (HDC) has constant probability density along its path and occupies the smallest possible volume in the variable space for a given probability content.We compute the contour using numerical integration based on a grid, i.e. we iteratively find the minimum probability density, f m , which leads to a contour containing the most likely environmental states which together have a probability of 1 − α.Defined this way a T -year environmental contour is exceeded on average every T years anywhere along the contour.This means precisely that such an environmental state is realized anywhere outside the environmental contour (and not in a further limited exceedance region).Highest density contours can be computed based on any probability density function, e.g. standard parametric sea state models, nonparametric models or extreme value models.The method's clear definition in terms of exceedance probability, α, as well as its straightforward computation makes it an attractive alternative to the established IFORM approach.

NomenclatureαFµ 2 , σ 2 mp
Exceedance probability [-] α W , β W , γ W Parameters of a Weibull distribution [-] f Cell-averaged joint probability density [-] F (f m ) Probability enclosed by a contour of f m probability density [-] fX1 Cell-averaged probability density in dimension 1 [-] fX2|X1 Cell-averaged probability density in dimension 2 conditional on x 1 [-] β Radius in U -space used in IFORM [-] Failure region [-] Parameters of a normal distribution [-] θ Angle [deg] µ Hs , σ Hs Parameters of a log-normal distribution [-] a 1 , a 2 , a 3 , b 1 , b 2 , b 3 Fitted parameters of the conditional distribution [-] C Set making up the environmental contour [-] Minimum probability density of the enclosed region / constant probability density along the contour [-] f * Normalized minimum probability density [-] H s Significant wave height, random variable [m] h s Significant wave height, realization [m] H s,25 25-year return value of the significant wave height based on its marginal distribution [m] H s Maximum significant wave height along the contour [m] a c c e p t e d v e r s i o n j Dimension index [-] K, L Number of grid cells in the respective dimension [-] k, l Grid cell index [-] LN () Log-normal distribution [-] M Random variable in general variable space [-] m Realization of the random variable in general variable space [-] n Total number of environmental states in a given time period [-] Number of variables / dimensions [-] P f Failure probability [-] p Hs Mixture coefficient [-] P r() Probability function [-] R Set enclosed by the environmental contour (highest density region) [-] Spectral peak period [s] T z Zero-upcrossing period, random variable [s] t z Zero-upcrossing period, realization [s] T z Maximum zero-upcrossing period along the contour [s]

Figure 1 :
Figure 1: Concept of an environmental contour.(a) The environmental contour encloses all variable combinations which must be considered in the design process (the design region).(b) Flowchart describing the design process utilizing an environmental contour.

Figure 2 :
Figure 2: Different definitions of environmental contours and their basis in univariate probability distributions.(a-c) Top: Univariate probability distributions (p = 1).Bottom: Example data and contours based on two-dimensional joint probability distributions (p = 2).(a) One-sided exceedance.(b) Two-sided exceedance.(c) Highest density regions (HDRs) with a minimum probability density, fm.

Figure 3 :
Figure 3: Established methods to calculate environmental contours.(a)Inverse first-order reliability method (IFORM)[34].The contour is defined as a circle in standard normal variables, U j .The points along the circle have to be transformed to the original environmental variables, X j .(b) Huseby et al.[17] Monte Carlo contour.The contour is computed with the environmental variables, X j , directly.(c) Constant probability density contour described by Det Norske Veritas[8].By its definition it is unclear how much probability is enclosed by the contour.(d) Jonathan et al.[21] constant exceedance probability contour.The calculation can be done in a general set of variables, M j , e.g. in the X-or U -space.In comparison to the other methods a different definition for the exceedance region is used (compare shaded areas).
a c c e p t e d v e r s i o n

Figure 4 :
Figure 4: Computation of the highest density contour (HDC) using a numerical grid.Shaded area = HDR, outline = HDC.(a)The variable space is discretized in equally sized grid cells and the average probability density, f , is calculated for each cell.The probability enclosed by a HDC of fm probability density is calculated by first finding all cells whose probability density, f , is greater than or equal the minimum probability density, fm, and then summing up the individual probabilities of these cells.(b) An environmental contour is computed by iteratively finding the minimum probability density, fm, that satisfies F (fm) = 1 − α.

Figure 5 :
Figure 5: Computed highest density contours.Along the contour probability density, fm, is constant and the enclosed region has a probability of 1 − α with α corresponding to a given T -year return period (T = 1, 10 or 25 years).Grid cell size is 0.05 m×0.05 s.

Figure 6 :
Figure 6: Expansion of the highest density contour.(a) The probability enclosed by the contour, F (fm), is 1 at a minimum probability density of fm = 0 and monotonically decreases to F (fm ≈ 0.12) = 0. Probabilities corresponding to the 1-,10-and 25-year contour are shown.The inlet illustrates the definition of F and fm.(b) Maximum variable values along the contour, H s and T z , as a function of minimum probability density, fm.The inlet illustrates that there is no (H s ,T z )-sea state along the contour.Instead, the (H s , Tz)-sea state has a Tz value different from T z and vice versa.

Figure 7 :
Figure7: Grid independence study.Quadratic grid cells with sizes ranging from 0.01 to 10 units of grid cell length are tested to evaluate grid convergence.(a) The contour's minimum probability density, fm, for a given return period is sensitive to grid cell size ∆Hs, ∆Tp.Sensitivity increases with grid cell length.(b) If grid cell size is too big minimum probability density, fm, can be half or double than the converged minimum probability density.Plotted is f * m which is minimum probability density, fm, normalized by the converged fm value.(c) Aiming for grid convergence with an error of less than 1 % we use grid cells with dimensions of 0.05 m×0.05 s (marked with a vertical line).

2 IFORMFigure 8 :
Figure8: Comparison of environmental contours derived using different methods.(a) Sketches showing expected differences in contour size due to different definitions.Some contours are defined in such a way that the maximum value along the contour, H s 25 , is equal the return value of the marginal distribution, H s25 , (middle).The highest density contour (HDC), however, is defined to enclose 1 − α probability and thus has a maximum value along the contour which is higher than the return value of the marginal distribution, H s 25 > H s25 , (right).(b) Sketch illustrating an IFORM contour and possible failure regions of a linear system of three components, F 1 , F 2 , F 3 .Since the contour contains less than 1 − α probability the system's failure probability can be greater than α.(c) A total of n = 73050 sea state data points have been Monte Carlo simulated representing a 25-year data set (scatter plot).The 25-year HDC (solid line) and the 25-year IFORM contour (short dashes) have similar shapes, but as expected the HDC is bigger.The 308.8-year IFORM contour or 25-year equi-shape contour (long dashes;[23]) encloses the same amount of probability as the 25-year HDC.(d) Comparison of maximum values along the contour, H s and T z .As expected by the different definitions, the HDC has the highest maximum significant wave height, H s , and maximum zero-upcrossing period, T z .

Figure 9 :
Figure 9: Environmental contours for mixture models.(a) Model 1 has a normal Tz-distribution, N (µ 2 = 10 s, σ 2 = 2 s), which smoothly blends into the original Tz-log-normal distribution.The highest density contours and IFORM contours have similar shapes.(b) Model 2 has a normal Tz-distribution, N (µ 2 = 15 s, σ 2 = 0.5 s), which leads to a second probability density maximum.Consequently, the highest density contour encloses two separated subregions.Due to its definition IFORM, however, encloses a single contiguous region.

mFigure A. 10 :
Figure A.10: Computer program to derive a highest density contour.Left: Code snippet written in the Matlab programming Right: Corresponding flowchart.