Analysing multivariate extreme conditions using environmental contours and accounting for serial dependence

Accurate statistical description of extreme environmental conditions is needed for risk assessment and management of marine structures and is a crucial input to design of any structure that need to withstand loads from environmental forces. Such descriptions are essentially multivariate extreme value problems, where the environmental loads are due to concurrent extreme combinations of several environmental variables. Typically, in coastal and ocean engineering applications, the simultaneous joint behaviour of significant wave height and wave period is of particular interest and is needed to describe the wave loads on marine structures. Environmental contours are often used to explore the extreme wave loads, and essentially consider extreme combinations of simultaneous significant wave height and wave period, and are based on a joint statistical distribution fitted to relevant metocean data. However, typical applications of environmental contours do not account for temporal dependencies in the environmental variables, and this may lead to an overestimation of extreme conditions. In this paper, an approach for partially accounting for serial dependence in the construction of environmental contours is proposed, based on simulating time-series of a primary variable which preserves both its marginal distribution and auto-correlation structure. It is shown that this gives lower estimates of extreme environmental conditions compared to conventional application of environmental contours that do not account for serial dependence. Hence, more accurate description of the extreme environment can be available for design and construction of structures exposed to environmental loads.


Introduction and background
Environmental contours are used to identify extreme environmental conditions governed by the simultaneous behaviour of several correlated environmental variables.Typical examples are extremes related to the joint behaviour of waves, wind and current for ocean and coastal engineering applications.Hence, environmental contours are tools for multivariate extreme value analysis often used in long-term extreme response analysis of marine structures, and represent a simple and approximate alternative to more computationally demanding full long-term analyses.One of the main advantages is that the structural response is de-coupled from the environmental description so that only a limited number of short term response calculations are required for a long-term analysis.That is, the environmental contours describe critical conditions for which short term analysis is performed [1].
Typically, environmental contours will be based on a set of environmental data, for example metocean data, and the application of environmental contours often consists of four main steps.The first step relates to collecting and preparing such data.This also include extreme event.In the bivariate case, a quantile will typically be vectorvalued and can be expressed on the following form [11]  , (, ) = { (, ) ∈ R 2 ∶   (, ) =  } , (1) where   (, ) corresponds to some exceedance or non-exceedance event.Different definitions of such events yields different types of extremes.Some examples are illustrated in the bivariate case in [12][13][14], i.e., cases where one of the variables is extreme, where both are extreme and where the value of one variable is extreme conditioned on the value of the other variable.Environmental contours represent other definitions of multivariate extremes, often described in terms of exceedance hyperplanes.However, there exist different definitions of environmental contours that do not refer to hyperplanes.Even with a given definition for multivariate extremes, corresponding to a specific choice of   (, ), there will not be a unique solution to the extreme quantile problem, and there is a continuum of solutions to the extreme value problem in the multivariate case.For example, in the bivariate case, if a solution ( 0 ,  0 ) to the extreme value problem exist, one could always find other solutions by slightly perturbing  0 and  0 so that Eq. ( 1) is still fulfilled [15].Thus, the solution will be a curve in the variable space associated with a non-exceedance probability defined in some way.
Environmental contours are one way of defining multivariate extreme sets.Environmental contours were initially proposed in [16,17] as equi-density lines, and the definition of environmental contours has continued to be an area of active research.Contours based on the inverse first order reliability method (IFORM) were introduced in [18,19], and remain the most commonly used approach.Such contours are based on exceedance hyperplanes in standard normal space and transformed to physical variable space.More recently, environmental contours defined in terms of exceedance hyperplanes in the physical variable space based on Monte Carlo sampling were proposed in [20,21].These contours, often referred to as direct sampling contours, have well-defined probabilistic properties, and will in some cases be quite different from the IFORM-contours, see the comparison studies in [22,23].In particular, the direct sampling contours will always be convex [24].Other approached for defining environmental contours exist [e.g., 15,[25][26][27][28][29][30][31], and several studies have investigated uncertainties associated with different aspects of environmental contours [23,[32][33][34][35].In this paper, direct sampling contours will be used, which correspond to exceedance regions in the form of exceedance hyperplanes.According to common industry practice, environmental contours are often based on a fitted joint distribution, fitted to all data and ignoring serial correlation [2,36].
Most applications of environmental contours are for twodimensional problems, but extensions to higher dimensions are in principle straightforward, see [37][38][39].However, the study presented in this paper is restricted to two dimensions and considers environmental contours for significant wave height (  ) and zero up-crossing wave period (  ).
Given that environmental contours have been an active area of research in recent years, a benchmark study was recently announced where practitioners and researchers of environmental contours were provided with metocean data and invited to construct contours for comparison [40].Several contributions were submitted for this exercise; results are presented in [41].One important observation from this benchmarking exercise is that serial correlation in the metocean variables is often neglected in the construction of environmental contours, and that this may introduce notable biases.This point was further elaborated in [42], where it was shown that neglecting serial dependence will lead to over-estimation of return values.Hence, this paper offers an approach to take serial correlation into account in the construction of environmental contours based on the direct sampling approach.This is done by establishing a statistical model that accounts for serial dependence in the primary variable, which will yield more accurate estimates of the extremes.
In univariate extreme value analysis, the effect of serial correlation is often accounted for by modelling the distribution of peak values, either using a block-maxima or peaks-over-threshold (POT) approach [43].For the latter method, de-clustering techniques must be applied to ensure that the samples from the tail are reasonably independent and identically distributed (IID) (see e.g.[43,44]).However, with limited amount of data, such approaches are wasteful since it only uses a very few extreme data points in the analysis and disregard most of the data.With the statistical model presented in this paper, all the data can be exploited and the effect of serial dependence is accounted for to reduce the bias.It is believed that this represents an advancement in state of the art, with the benefit of more accurate extreme value estimates.
Identifying peak events is not straightforward in the multivariate case, where combinations of the variables that may not be extreme in any one variable may be regarded as extreme when their joint probability is considered.Indeed, environmental contours typically range over both extreme and non-extreme values of all variables involved (including extremely small values).Moreover, the interest is in environmental loads and responses, which means that the concurrent behaviour of the environmental variables is needed.Hence, some multivariate extreme value methods based on component-wise extremes will not be applicable [45,46].Some methods, such as the bivariate ACER model [47], explicitly account for the serial correlation in bivariate data.
An approach for defining environmental contours that accounts for serial correlations in the observations was recently proposed in [48].The method does not require fitting a joint distribution model.Since defining exceedance regions in terms of hyperplanes is equivalent to defining univariate exceedance regions under rotations of the axes [49], the problem of estimating environmental contours can be reduced to a series of univariate problems under various rotations of the axes.A standard univariate POT approach with de-clustering is then applied to the rotated data, to account for serial correlations.
This paper proposes an alternative approach to take serial correlation into account in the construction of environmental contours, based on a direct sampling approach.This approach to construct environmental contours is based on Monte Carlo samples from the joint distribution of environmental variables, where the data is projected in selected directions and hyperplanes are defined in all directions based on the target exceedance probability.The environmental contours are then constructed based on the intersection points of these hyperplanes; see [20,21] for details.Hence, this paper proposes to construct environmental contours on simulated data from a time-series model that preserves the marginal distribution and the serial correlation in the data.A conditional modelling approach, in line with current industry practice, will be employed where the zero up-crossing wave period is modelled conditionally on significant wave height, but significant wave height will be modelled by a time-series model with a fitted marginal distribution and temporal correlation rather than independently from the marginal distribution.In this way, environmental contours that accounts for serial correlation of significant wave height can be constructed.The main challenge with this approach is to establish a statistical model that allows simulation from such a time series with the desired properties.
The remainder of the paper is organized as follows: Section 2 reviews how return periods and return values are defined, and discusses the biases caused by neglecting serial correlation.Section 3 presents the metocean data used in this study and Section 4 presents the statistical models fitted to the data.This includes the time-series model with desired marginal and correlation structure for   as well as the conditional model for   .Section 5 presents the resulting contours, and a discussion of the results are given in 6.Finally, Section 7 gives a summary and some concluding remarks.

Return periods and return values for serially correlated data
Informally speaking, the return period of a value of a parameter is the average time between exceedances of that value.The value corresponding to a specific return period is called the return value.To understand the effect that serial correlation has on estimates of return periods and return values, consider two variables  and X, with the same marginal distribution,   , but where observations of  have some level of serial correlation and observations of X are independent.Serial correlation effectively reduces the number of independent observations in a sequence of a given length, and therefore reduces the probability of large values.Since both  and X have the same marginal distribution, but there is a lower probability of large values in  in a given period, we see that defining return values in terms of the marginal distribution of all observations,   , cannot be correct when serial correlation is present.
For serially correlated data return periods can be defined in terms of the distribution of peak values in the time series [43].Suppose that local peaks,   , have been identified in a time series, using some declustering criteria, such that consecutive observations of   are approximately independent.Then the return period of level  is where    is the distribution of the peaks and  is the average number of peaks per year.This definition is also applicable when annual maxima are used instead of local peaks.In this case    is replaced by the distribution of annual maxima and we set  = 1.Alternatively, the  -year return value can be defined as the solution of [50] where   is the distribution of the maximum value of  in a  -year period.It can be shown that definitions (2) and (3) are asymptotically equivalent (see e.g.[42,50]).For the sequence X, since all observations are independent, the return value of level  can be defined using (2) as where  is the (fixed) number of observations per year.It was shown in [43] that the ratio of return periods for the independent and dependent sequences is less than or equal to 1.That is The ratio   quantifies the bias introduced when ( 4) is used to define return values for serially correlated data.It can be interpreted as a sub-asymptotic extremal index, which converges to the usual extremal index, which by definition is an asymptotic quantity, used to quantify the effect of serial correlation in extreme events (see [43,51] for details).

Metocean data description
The data used in this study are one of the datasets provided for the recent benchmarking study on environmental contours, i.e., dataset A, from [40].Several environmental contours established based on these data are presented in [41].The data contains hourly observations of significant wave height (  ) (in ) and zero upcrossing wave period (  ) (in ) over a 10-year period from 1996 to 2005.It is noted that 10 years of data is rather short for robust extreme value estimation, and uncertainties will be large see e.g. the discussions in [52,53].
A scatterplot, trace plots and empirical densities for these data are shown in Fig. 1, along with summary statistics in Table 1.These plots and summary statistics reveal that   and   are positively correlated, The traceplots illustrate that there are notable serial dependencies in the data.For further details and description of the data, reference is made to [40,41].For reference, results from two crude marginal extreme value analyses are shown in Table 2, where the 10-year return value is estimated for both variables.The first is obtained by fitting a generalized Pareto distribution to peaks over threshold data with a threshold of, respectively, 5 m for   and 11 s for   and a cluster separation distance of three days (denoted as GPD in the table).The second is based on fitting a generalized extreme value distribution to semi-annual maxima (denoted as GEV in the table).

Statistical modelling
For the purpose of structural reliability analysis and design of marine structures, the joint distribution of governing environmental variables such as significant wave height and wave period is often modelled by a conditional model.Thus, one variable is considered the primary variable, typically the one with the largest influence on the environmental load.A joint distribution is constructed based on fitting a marginal distribution to the primary variable, and then conditional models to the other relevant variables [3][4][5].Such a model was fitted to the data used in this study, as reported in [54], but without considering serial correlation.
In order to be able to construct environmental contours that account for the serial correlation in the data using the direct sampling approach, one may fit a time series model with appropriate marginal distribution and temporal correlation structure.This is not straightforward, and there are different approaches that could be used [55].For example, [56] proposes a diffusion-type model with given marginal distribution and autocorrelation function, and [57] describes a diffusion model with Weibull marginals for modelling wind speed.Models based on stochastic differential equations are presented in [58,59].Weibull and gamma autoregressive processes are presented in [60], and Gaussian and non-Gaussian autoregressive models are discussed in [61].However, in this paper, an approach based on transforming a Gaussian parent autoregressive time series to have the desired marginal distribution and temporal correlation structure will be taken, as suggested in [62]; see also [63][64][65].In the following, a recipe for constructed such a model is outlined.A joint statistical model is constructed by fitting a time-series model that preserves the marginal distribution and the correlation structure for the primary variablesignificant wave height -and a conditional model for zero up-crossing wave period conditioned on significant wave height.It should be noted that with this approach, the serial correlation of   is not modelled directly but is a result of a conditional distribution conditioned on   and would only be partially preserved.
The time-series of interest, (), has a marginal distribution,   () and an autocorrelation function,   () and it is assumed that this has a parent Gaussian time-series, () with a standard Gaussian marginal   () and another autocorrelation function   ().Hence, modelling and simulation of the time-series () correspond to modelling and simulating the parent Gaussian time-series and finding the transformations that transform the parent time-series to the desired one, i.e., determine functions  and  such that Simulating from a Gaussian time-series is straightforward, and the normal distribution have several well-known properties (see e.g.[62,66]).Hence, if one can find a parent Gaussian time-series that can be transformed into the time-series (), simulation of () is straightforward.As pointed out in [62],  is easily identified as ) , which will give the desired marginal distribution.However, this transformation of the marginal distribution does not preserve the autocorrelation structure, which depends on the marginal distribution   ().Hence, one needs to establish the autocorrelation structure of the parent Gaussian time-series,   () that yields the required autocorrelation structure for the target time-series () after transformation.
In the bivariate case, it can be shown that for an arbitrary transformation, the correlation between the transformed variables will be smaller than the correlation between the initial standard normal variables (see [62]), i.e., meaning that the   values needs to be inflated to obtain the target   .Moreover, the correlation coefficient of two lagged variables ( 1 ) and ( 2 ) can be expressed as where  is the expectation operator,  =  2 −  1 and   and   are the mean and the standard deviation of .Now, assuming that where (, ; ) denotes the density function of the bivariate standard normal distribution with correlation .This gives a relationship between   () and   () and can be used to calculate   () from known   (), i.e.
Note however, that in order to model the time series with known   (), one needs the inverse transformation to find   () from   ().Even though the double integral above does not have an analytical solution in general, its numerical estimation is straightforward.Hence, the modelling problem is solved by defining a correlation transformation function to estimate the autocorrelation of the parent Gaussian time-series from a given target autocorrelation structure that may be estimated from the data.The following parametric form of such an autocorrelation transformation function is proposed in [62]: This can be used to estimate a parametric autocorrelation function of the parent Gaussian time-series,   () based on a parametric autocorrelation function for the target time-series,   ().The parameters  and  can be estimated by calculating   for   = (0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95) according to Eq. ( 10) and then conducting a least-squares minimization of the assumed model (11) to the set of (  ,   )-points (see [62] for details).With this approach, modelling a time series with desired marginal and correlation structure involves fitting a marginal parametric distribution and estimating the empirical autocorrelation structure from the time-series, fitting the parametric autocorrelation transformation function to obtain the autocorrelation structure of the parent Gaussian time-series, fit an autoregressive model, for example by solving the Yule Walker equations [e.g., 67], to generate the Gaussian parent time-series and transforming the Gaussian time-series using the transformation Hence, in the following, environmental contours for significant wave height (  ) and zero up-crossing wave period (  ) will be based on simulated time-series for   with prescribed marginals and autocorrelation structure and concurrent values for   simulated as independent values from a conditional model for   given   .Therefore, any serial correlation in   is assumed to be, at least partly, accounted for by the serial correlation in   and the dependence structure between   and   .The validity of this assumption is discussed in Section 6.

Joint distribution
First, a joint distribution will be fitted to the data disregarding the serial correlation.In this study, a conditional model with a threeparameter Weibull distribution for   and a conditional log-normal model for concurrent   will be assumed.That is, where and This is the same model that was used in [54], where environmental contours were established for this dataset without considering serial correlation.Fitted model parameters are presented in Table 3.It is noted that this particular model was fitted to the data by way of minimizing the second order Anderson-Darling statistic.This is emphasizing on the upper tail, and may not fit the body of the distribution as well as more well-known methods as the maximum likelihood or the method of moments (see e.g.[44]).

Temporal correlation structure
The temporal correlation structure is estimated from the empirical time series.Then, a number of alternative parametric autocorrelation functions are fitted to this, and the one providing the best fit is selected for further analysis.There are many parametric autocorrelation functions to choose from, see e.g.[62] and in this study, six different parametric functions are tried out, i.e. the Weibull, the Burr, the logarithmic, the fGn, the generalized fGn (gfGn) and the Pareto autocorrelation functions.The empirical autocorrelation function is shown together with the fitted parametric candidate models in Fig. 2. The functions are fitted by least squares, and the residual sum of squares (RSS) are indicated for each candidate model in the plots.It is observed that the 3-parameter Burr-type autocorrelation function fits best to the data, and this would typically be preferred.However, it is noted that this autocorrelation function is known to sometimes yield not positive definite covariance matrices (see [68]), and therefore this model may not always be appropriate.The preferred autocorrelation function must be determined on a case-by-case manner, according to how well it fits the data, and it is also possible to use the non-parametric empirical autocorrelation function.The parametric forms of the autocorrelation functions explored in this study are given in Eq. ( 16).

𝜌 𝐵𝑢𝑟𝑟 (𝜏; 𝜁 , 𝜂, 𝜃)
It can be observed that most of the parametric autocorrelation functions fit reasonably well to the data for lags up to 100, except for two of the parametric alternatives.It is also noted that the choice of autocorrelation function may influence return value estimates, as seen in [65], so some efforts should be made to find the best possible fit to the data.However, if a good fit cannot be found, it will be possible to simply assume the non-parametric autocorrelation observed in the data.

Simulated time-series
Having estimated all model parameters, one may simulate timeseries from the statistical model corresponding to any desired number of years.The simulated data would then have the desired marginal distribution and the desired autocorrelation structure.One may also simulate from the marginal distribution assuming IID.Fig. 3 shows about 1 year of empirical data for   and compares to simulated data based on the time-series model and simulated data assuming IID.The bottom plot in Fig. 3 contains first IID simulated data corresponding to half the empirical data duration and then simulated time-series data for the remainder of the empirical data length.It is clearly seen that the data simulated when accounting for the serial correlation better resemble the observed data.
One may also compare the densities of simulated data with the density of the empirical data, as shown in Fig. 4, where 500 time series of 10 years' length are simulated, assuming IID and with the time-series model that preserves the autocorrelation, respectively, and the simulated densities are compared to the empirical density of the observed data.The densities have been fitted based on a goodness-offit measure that focuses on the upper tails of the distribution, i.e., the second order Anderson-Darling statistic [44,69], and therefore there is a notable mismatch for small values of   .However, since the upper tail of the significant wave height distribution is more important for structural reliability of marine structures, the poor fit on the lower tail is acceptable.It is noted, however, that better fit of this part of the distribution could presumably be obtained if other fitting methods were employed, and this is recommended in applications where the lower tail of the distribution is of importance.Exceedance probabilities (on logarithmic scale) are also compared to highlight tail behaviour.In the upper plots, the densities and exceedance probabilities are estimated from the complete simulated datasets of 5000 years, whereas the lower plots show the densities and exceedance probabilities for individual 10-year periods of simulated data.The top plots clearly demonstrate that the marginal distributions are the same, regardless of whether data are simulated IID from the fitted distribution or from the timeseries model with the same marginal distribution.However, there are some differences in the simulated individual 10-year datasets, due to sampling effects.The sampling uncertainty associated with the largest observation in a sample is very large.As the observations in a sample are random quantities, the return periods associated with each observation are also random quantities.It can be shown that a 1 − 2 confidence interval (CI) for the return period associated with the largest observation in an -year time series tends to (−∕ log(), −∕ log(1 − )) as  → ∞ [70], and the approximation is still reasonable for  = 10 years.So, for individual 10-year time series, a 95% confidence interval for the return period associated with the largest observation would be (2.7,395) years.The width of the confidence interval decreases as the ratio ∕ increases, where  is the length of the time series and  is the target return period (see [70] for details).Fig. 5 shows the distribution of the 10-year maxima from the 500 simulated 10-year time-series.The 10-year return values estimated from these distributions are also shown, where the return values are estimated as the (1 − 1∕10) 10 ≈ 0.3487 quantiles of the maximumvalue distributions according to Eq. (3).It is observed that the return values for the 10-year maxima are higher for the IID data compared to the serially-correlated time series.The vertical lines for the combined data correspond to ignoring serial correlation and estimating the return value according to Eq. ( 4) as the quantile from the combined simulated datasets from the IID and TS simulations, respectively.It is observed that these generally agrees with the return value estimate from the distribution of 10-year IID datasets, apart from small effects due to sampling variability, and the difference can hardly be recognized in the figure .It is also interesting to compare the return value estimates from the simulated data with the summary statistics and the estimated return values from classical extreme value methods presented in Tables 1-2: 10-year return value estimated from simulations assuming independent and identically distributed data is 8.967 m whereas 10-year return value estimate from simulated time series is 7.229 m.It is observed that the return value estimated based on the simulated data from the time series model is much closer to the return value estimates from the classical extreme value analysis, and also closer to the maximum value observed in the data.Hence, the positive bias from ignoring serial correlation is avoided.4); the other return value estimates are based on the distribution of 10-year maxima, for simulated data assuming IID and the time-series model, respectively, according to Eq. ( 3).The maximum observed values are also indicated by vertical line segments, illustrating that accounting for serial correlation effectively corresponds to simulating shorter time periods compared to IID.

Environmental contours accounting for serial correlation
Environmental contours are often used for describing extreme environmental conditions.With current industry practice, such contours are typically constructed from exceedance probabilities based on Eq. ( 4), ignoring the effect of serial correlation.An alternative approach to construct environmental contours that takes the effect of serial correlation into account is proposed.With this approach, a total of   -year contour estimates can be obtained by simulating   -year timeseries from the statistical model and estimate the contours for the most extreme conditions in each individual simulated  -year period of data.This will give a Monte Carlo distribution of  -year contours which can be used to construct environmental contours corresponding to year return values according to Eq. ( 3).The main steps involved in estimating  -year environmental contours using this approach are as follows assuming an appropriate time-series model has been established (see [20,21] for further details on how direct sampling contours are constructed): 1. Simulate   -year time series from the statistical model 2. For selected angles of rotation , get the  angular maxima from the simulated time-series to estimate the distribution of  -year angular maxima,   (  ), where   denotes the data projected in direction :   =  cos  +  sin  for bivariate data (,  ), see, e.g., [20,21].3. Calculate return values, in each direction, according to Eq. (3) 4. Calculate  -year environmental contours from the angular return values by the direct sampling approach [20,21] The results of simulating  = 500 10-year time-series from the statistical models described in Section 4 and calculating environmental contours for the maxima in each individual time-series are shown in Fig. 6. Results are shown for both data simulated from the joint model assuming IID and for data simulated from the time-series model accounting for autocorrelation.The results from the standard way of constructing environmental contours (ignoring serial correlation) based on the combined 5000-year datasets are also indicated in these plots (referred to as ''combined'').Contours corresponding to the 10-year return values, calculated from Eq. (3) based on the simulated 10-year maxima, are also included in the figures, for both the IID and the timeseries simulated data, together with 95% confidence bounds of the 10-year maxima.Vertical lines indicate the corresponding values for the estimated 10-year return values for   by the different approaches.It is observed that the results from the standard approach on the combined data are the same as the results for the 10-year return value based on individual 10-year datasets assuming IID.However, a notable difference is observed for contours based on time-series simulations taking autocorrelation into account.This corresponds to the bias in neglecting serial correlation, and these results indicate that this bias is notable and may have practical consequences in design and assessment of marine structures.
The environmental contours corresponding to the 10-year return value are compared in the same plot as shown in Fig. 7.The contours are based on the individual simulated 10-year periods and calculated according to Eq. ( 3), assuming IID and autocorrelation.Again, a notable difference can be seen, and in particular, contours accounting for serial dependence correct for the positive bias under the IID assumption.This is most pronounced in the primary variable, significant wave height, and not in the secondary variable   .This is presumably because the serial correlation in the secondary variable is not well captured by the statistical model (see also the discussion in Section 6 below).Notwithstanding, it is suggested that the approach outlined in this paper can be adopted when calculating environmental contours as one way of accounting for serial dependence in observed data.
An alternative way of estimating environmental contours based on time-series simulations would be to, rather than simulating  year time series and calculating  -year contours from a distribution of simulated  -year maxima according to Eq. (3), to simulate  oneyear time series from the statistical model and then calculate  -year contours from the distribution of simulated annual maxima according to Eq. ( 2).Typically,  >  but since the individual time-series would be much shorter this could be computationally more efficient.Moreover, the simulated data could be used to estimate environmental contours for any return period  > 1.However, the overall idea is the same and the results would presumably be similar.The uncertainty of the return periods associated with the contours could then be calculated explicitly as a function of  and  as described in [70].Fig. 8 show the distribution of maxima from  = 1000 simulated 1-year time series, assuming IID data and accounting for the autocorrelation, respectively.Vertical bars in the plot indicate the estimated 1-year return values from Eq. ( 3) as well as estimated 10-and 100-year return values from Eq. ( 2).It is clearly seen that the IID assumption leads to a positive bias in the return value estimates.The corresponding 1-, 10-and 100-year environmental contours based on the  = 1000 simulated yearly time series are shown in Fig. 9 One advantage of the proposed compared to the approach proposed in [48] is that no de-clustering is required.However, a disadvantage is the need for fitting a time-series model, and the results will be depend on how well such a time-series model capture the temporal correlations.An alternative approach would be to use a block-resampling approach, rather than a time-series model, as described in [71].The same method for constructing contours as described above could be used for time series generated from any model or resampling approach.

Discussion
This paper has presented a way of defining environmental contours that take serial correlation properly into account.Essentially, time Fig. 6.Contours for the 10-year return value based on simulated time series; assuming IID (left) and accounting for autocorrelation (right).The plots show M = 500 individual environmental contours calculated from independent 10-year simulations as well as contours corresponding to the 10-year return value from all the simulations (dashed lines).Direction-wise 95% confidence intervals are included (dotted lines).Contours calculated from the combined datasets ignoring serial correlation are also shown in the Figures.Vertical lines correspond to return value estimates of the primary variable,   and show results from the combined dataset (solid lines) as well as results based on the distribution of 10-year maxima (dashed lines) from simulations assuming IID and accounting for serial dependence, respectively.series are simulated by transformation of a parent Gaussian autoregressive model to yield the target marginal distribution and serial correlation.It has been demonstrated that this yields different contour estimates compared to the commonly used IID assumption, and that the positive bias in IID contours can be avoided.Even though this is believed to be an improvement, some aspects have not been accounted for.
In the statistical models used to construct the environmental contours, seasonal effects have been ignored.Possibly, such effects could also be accounted for in more advanced time-series models, e.g.assuming seasonal ARIMA models rather than simple autoregressive models as the parent time-series.Alternatively, seasonal effects could be taken into account by preprocessing of the data before the statistical models are estimated, as suggested in e.g.[13].Notwithstanding, seasonality  has not been considered in the current work, and remains an area for further research and potential for improvement of the modelling.Other deterministic components of the observed time-series, for example possible long-term trends due to climate change [72,73], should also ideally be filtered out before the autocorrelation function is estimated [74][75][76] but this is left for further work.If climate change can be expected to change the future environment, one may model this explicitly in the statistical models, or one may want to add additional safety factors accounting for this to ensure the necessary conservatism and safety level of the design [77].At any rate, the effect of climate change can potentially be included in the statistical models, for example as an additive deterministic component.
It is noted that only 10 years of data have been available for this study, and that this is generally too short to give robust estimates of extreme conditions; there will be a need to extrapolate the probability distributions beyond the support of the data.Hence, the actual 10-year return values is not known, making it difficult to evaluate directly the impact of accounting for serial correlation on the extreme estimates.Notwithstanding, it is known that neglecting serial correlation yields a positive bias, and it is observed that return value estimates without accounting for serial dependence gives higher return value estimates than the ones obtained from the time series model accounting for it.Therefore, it is assumed that accounting for serial dependence gives an improved extreme value estimate that is closer to the true value, and that corrects for this positive bias.However, exact quantification of this improvement is difficult without longer datasets and knowledge of the true return values.
The statistical model assumed in this paper consists of a timeseries model for the primary variable,   , which preserves the desired marginal distribution and autocorrelation structure, and a conditional distribution for the concurrent value of   .It is assumed that this is sufficient and that preserving the autocorrelation structure of the primary variable as well as the dependence between the primary and secondary variable by way of the conditional distribution yields improved description of the extreme environmental conditions.However, it is acknowledged that the autocorrelation structure of the secondary variable will not necessarily be preserved.Hence, further research can be directed towards time series models that preserves all marginals,  the joint distribution and the autocorrelation structures of multivariate time-series (see, e.g., [78][79][80]).The empirical autocorrelations are shown for the empirical and simulated data in Fig. 10 for both the primary and the secondary variable.It is observed that the autocorrelation is preserved reasonably well for the primary variable, but not for the secondary.It is observed that some serial correlation is obtained from the conditional modelling approach, but that this is far less than in the observed time-series.Notwithstanding, in practical engineering applications, the extremes of the primary variable, significant wave height, will be most important, so it is believed that the model proposed in this paper represent an improvement that is practically significant in ocean engineering applications.
When the interest is in the extremes, it could be argued that it is not necessarily the standard autocorrelation function that is of most relevance, but rather the level-dependent serial correlation expressed by the conditional probability (, ) =   (( + ) > |() > ), which may be illustrated by extremograms [81,82].The correlation structure for high levels, , should also be preserved in the time series to ensure the effect of serial correlation on extremes are properly accounted for.In this study, the emphasis has been on simulating time-series that preserves the marginal distribution and the autocorrelation structure, rather than preserving the extremogram, and further studies are recommended to construct models that preserve also extremal serial correlation.However, it is possible to compare the sample extremogram of the data and the simulated time-series, as shown in Fig. 11 for simulated significant wave height for levels corresponding to quantiles  = 0.99, 0.999 and 0.9995.As can be seen in these plots, the simulated time series preserves the extremal correlation quite well, and this is clearly an improvement compared to assuming IID data.

Summary and conclusions
This paper considers the construction of environmental contours for the description of multivariate extreme environmental condition, a practice that is often used in design and assessment of marine structures.However, the effect of serial dependence is often neglected in engineering practice, and this paper proposes a way to construct contours while accounting for the serial dependence.A statistical model constructed by a joint conditional model for the variables and a timeseries model preserving the autocorrelation structure for the primary variable is established and applied to data of significant wave height and zero up-crossing period.This latter time-series model is constructed by transforming a parent Gaussian autoregressive model using appropriate transformations for the marginals and the autocorrelation structure.It has been demonstrated that environmental contours may be constructed based on such statistical models, using the direct sampling approach.Moreover, results indicate that the resulting environmental contours avoid the positive bias one normally gets when ignoring serial dependence, as is current industry practice.This again, means that a more accurate description of the extreme environment is achieved that can be used for more informed decisions in design and operation of marine structures.
Further research could investigate how to account for seasonal effects and the serial correlation of secondary variables as well as models that better capture the serial correlation of extremes.Notwithstanding, the proposed approach represent a significant improvement compared to current engineering practice of ignoring serial correlation and gives more accurate environmental contours.

Fig. 3 .
Fig. 3. Empirical (top) and simulated (bottom) time-series of significant wave height.Index represents time in hours.

Fig. 4 .
Fig. 4. Densities and exceedance probabilities of empirical and simulated data; based on combining all simulated data (top) and for individual simulated 10-year periods (bottom).

Fig. 5 .
Fig. 5. Distribution of extremes from individual 10-year sets of simulated data; IID and time-series simulations.The vertical lines correspond to different estimates of the 10-year return value: Combined data correspond to return value estimates based on the combined datasets of 500 × 10 years of data and ignoring serial correlation according to Eq. (4); the other return value estimates are based on the distribution of 10-year maxima, for simulated data assuming IID and the time-series model, respectively, according to Eq. (3).The maximum observed values are also indicated by vertical line segments, illustrating that accounting for serial correlation effectively corresponds to simulating shorter time periods compared to IID.

Fig. 7 .
Fig. 7. Comparing contours based on time-series modelling, assuming IID or accounting for autocorrelation in the data; based on simulated 10-year time series.

Fig. 8 .
Fig. 8. Distribution of maxima from simulated one-year time series.

Fig. 10 .
Fig. 10.Autocorrelation structure of observed and simulated time series for both primary and secondary variable.

Fig. 11 .
Fig. 11.Comparing the extremogram of the data with that of simulated time-series data for selected quantile levels of significant wave height.

Table 1
Summary statistics for the metocean data.

Table 2
10 year return value estimates from marginal extreme value analyses.
and that both marginal distributions are right skewed and leptokurtic, i.e., that they have fatter tails than a normal distribution.The 0.99989 quantile of the data would correspond to the one-year return value if not considering serial dependence, i.e. 1 − 1∕(365.24* 24) ≈ 0.99989.

Table 3
Fitted model parameters for the joint distribution of   and   .