Correlation of Concurrent Extreme Metocean Hazards Considering Seasonality

: Simultaneous occurrence of metocean variables can present a multihazard to maritime systems. However, simpliﬁed design approaches to assess simultaneous signiﬁcant wave heights and wind velocities are lacking, especially if seasonality is considered. This is addressed in this study by using extreme signiﬁcant wave heights and companion wind velocities recorded in the Gulf of Mexico. Time-dependent, generalized extreme value (GEV) models and classical regression are the basis to propose a simpliﬁed approach to estimate correlated extreme signiﬁcant wave heights and wind velocities associated with given return periods, accounting for seasonality and including measures of uncertainty. It is found that the proposed approach is a new but simple method to adequately characterize the concurrent extreme metocean variables and their uncertainty. It is concluded that the method is an e ﬀ ective probabilistic design tool to determine simultaneous extreme signiﬁcant wave heights and companion wind velocities for desired return periods and seasonality.


Introduction
Maritime structures, such as offshore facilities, breakwaters and other systems, are subjected to the effect of metocean variables, including extreme waves and extreme winds. These metocean variables can be considered as hazards that impose demands to these systems. Metocean variables can be characterized with a probability density function (PDF) frequently adopted from extreme value theory [1]. Among the metocean variables, extreme waves and winds can impact offshore facilities [2]; correlation of these two variables is important for design because they can occur simultaneously [3].
Recent studies on extreme waves and winds include the assessment of significant wave heights in the South China Sea using hindcast data [4], wind models to simulate waves in the South China Sea and the East China Sea [5], the use of radar images to estimate significant wave height [6], wave height forecasting during extreme events [7], assessment of wave heights by using generalized Pareto and Gumbel distributions [8], uncertainty assessment in extreme significant wave height by using different techniques [9] and a review on available wind-wave models and their limitations [10]. Additionally, hindcast for wind-wave for several decades was calculated and compared with observed data for estimating wave climate change and

Extreme Value Analysis Considering Seasonality
In this section the method to perform the extreme value analysis considering seasonality is briefly summarized, since details are described in a previous study [28]. In the next section, the model is to be first applied separately to the significant wave height and wind velocities.
The model employs the generalized extreme value (GEV) family of distributions [29]. Instead of characterizing a random variable with the Gumbel, Weibull or Frechet PDF or cumulative functions (CDF), a single mathematical expression can be used, since it can represent the three. This is the GEV model used here and given by: (1) or given by: if the PDF is to be used. In Equations (1) and (2), [a]+ implies the max (a,0), −∞ < µ < ∞ is the location parameter, ψ > 0 is the scale parameter and ξ the shape parameter; ξ leads to the Gumbel distribution when equal to zero, and to the Frechet and Weibull distributions when ξ > 0 and ξ < 0, respectively. For the GEV models of Hs and V w , the block maxima concept is used. The monthly maximum Hs values are considered as the samples for the extreme value analysis, while for the samples of V w , they do not necessarily correspond to the maximum, but to the companion wind velocities recorded when Hs occurred. The seasonality is a cyclic behavior linked to climatic patterns, which behave in a sinusoidal way every year. The common stationary model for extreme value analysis is extended to include time dependency. It is assumed that monthly maxima and companion values are independent random variables. Furthermore, it is assumed that the maximum (or companion) monthly values Z t of observed Hs and V w in month t follows the GEV distribution referred to above, which location µ(t) > 0, scale ψ(t) > 0 and shape ξ (t) parameters are a function of time.
The selective search methodology known as stepwise approach is used to find the best model, which encompasses the forward selection and backward elimination methods, which are described elsewhere [30]. The uncertainty of the selected model is accounted for by standard errors and confidence intervals [30]. The present study proposes an approach to compute concurrent metocean variables (Hs and V w ) associated with given return periods, for which the individual values associated with each variable are required. The values associated with return periods for each variable considered separately can be obtained with the following equation: X q (t, θ) = X q (µ(t), ψ(t), ξ(t)) = where q corresponds to the exceedance probability as defined by G t (x) = 1 − q and the estimated quantile, and X q (t,θ) corresponds to the time-dependent value associated with the return period T -yr = 1/q for a metocean variable of interest. The model described in this section is applied to Hs and V w in the following sections; first, for each individual variable, and then subsequently, the approach to correlate both metocean variables is presented.

Nonstationary Models Considered Separately
The time-dependent GEV model described in the previous section is first applied separately to each of the metocean variables. The significant wave height Hs and wind velocity V w described before are employed. Recall that while for the GEV analysis of Hs the monthly maxima from buoy 42,001 are used, the V w data do not necessarily correspond to the monthly maximum recorded values, but to the companion wind velocities that occurred simultaneously with the monthly maximum Hs. Figure 1 shows the maximum Hs recorded monthly (depicted with dots) and the 30 yr return period values, T-30 , estimated with the GEV model (depicted with a solid line). As in a previous study [28], clear peaks approximately around the seasons of December-February and August-October, associated with cold fronts and hurricanes, respectively, are observed in Figure 1.

Nonstationary Models Considered Separately
The time-dependent GEV model described in the previous section is first applied separately to each of the metocean variables. The significant wave height Hs and wind velocity Vw described before are employed. Recall that while for the GEV analysis of Hs the monthly maxima from buoy 42,001 are used, the Vw data do not necessarily correspond to the monthly maximum recorded values, but to the companion wind velocities that occurred simultaneously with the monthly maximum Hs. Figure 1 shows the maximum Hs recorded monthly (depicted with dots) and the 30 yr return period values, T-30, estimated with the GEV model (depicted with a solid line). As in a previous study [28], clear peaks approximately around the seasons of December-February and August-October, associated with cold fronts and hurricanes, respectively, are observed in Figure 1. To inspect whether a difference using the updated information in the present study (i.e., with data up to 2019) and that from the previous study [28] (with data up to 2012) exists, the T-30 values reported in [28] are also depicted in Figure 1 (dashed lines). It can be observed that inclusion of recent data points in the analysis does have an effect in the projected values, leading to a decrease in the peak Hs values for hurricane and cold front seasons-more noticeable for the former-and an increase at midyear months.
To further investigate the impact of selecting different time windows in the assessment of projected values associated with a given return period, two more lines are shown in Figure 1, which were computed using the oldest (dotted line) and the newest (dashed-dotted line) maximum 20 available recorded monthly data per year; selection of at least a 20 years range for the analysis is recommended to properly characterize metocean variables probabilistically [3]. There is a noticeable impact of using these different time windows to assess the Hs values for T- 30. Determining if this could be attributed to climate change is feasible, since studies in the literature [31] report that higher or lower metocean variables could be attributed to climate change, which depends on region and location, although it is important to acknowledge uncertainty in claiming any To inspect whether a difference using the updated information in the present study (i.e., with data up to 2019) and that from the previous study [28] (with data up to 2012) exists, the T-30 values reported in [28] are also depicted in Figure 1 (dashed lines). It can be observed that inclusion of recent data points in the analysis does have an effect in the projected values, leading to a decrease in the peak Hs values for hurricane and cold front seasons-more noticeable for the former-and an increase at midyear months.
To further investigate the impact of selecting different time windows in the assessment of projected values associated with a given return period, two more lines are shown in Figure 1, which were computed using the oldest (dotted line) and the newest (dashed-dotted line) maximum 20 available recorded monthly data per year; selection of at least a 20 years range for the analysis is recommended to properly characterize metocean variables probabilistically [3]. There is a noticeable impact of using these different time windows to assess the Hs values for T -30 .
Determining if this could be attributed to climate change is feasible, since studies in the literature [31] report that higher or lower metocean variables could be attributed to climate change, which depends on region and location, although it is important to acknowledge uncertainty in claiming any categorical trend [31]. This potential variation in wave climate could importantly impact maritime structures [31]. Future studies are desirable to further investigate this topic; nevertheless, it can be concluded that selection of different time ranges does have an impact in the projected return values of Hs (and other metocean variables). Moreover, it is accepted that using more metocean data leads to more reliable extreme value analyses [2,3], but the use of less data and more recent data could lead to a higher projected return period value for some months. Consequently, future research and discussion on this issue is recommended.  Figure 1, except that V w values are shown. Less accentuated peaks can be observed for hurricane and cold front seasons in Figure 2, perhaps because not the maximum, but the companion values are employed, even though the trends are similar, at least for the updated (solid line) and previously used [28] time ranges (dashed lines).  Figure 1, except that Vw values are shown. Less accentuated peaks can be observed for hurricane and cold front seasons in Figure 2, perhaps because not the maximum, but the companion values are employed, even though the trends are similar, at least for the updated (solid line) and previously used [28] time ranges (dashed lines).
For the scope of this study, only models corresponding to the solid lines in Figures 1 and 2 are used (i.e., the whole set of available data at buoy 42,001 is considered), although projected metocean variables for other return periods are also used. The values for other return periods, used later in this study, for each month are listed in Appendix A in Tables A1 and A2 for Hs and Vw, respectively. No details in this study are given about the parameters of the GEV models as were provided in previous studies (e.g., [28]), since attention is focused on proposing a simplified approach to correlate extreme Hs data and their companion Vw data associated with given T-yr values. This is developed in the following section.

Simplified Approach for Return Period Correlated Values
The proposal to estimate the simultaneous occurrence of extreme metocean variables is based on a modified version of the classical linear regression technique. This consists in including estimated values from the time-dependent GEV models in the previous section-together with the simultaneously observed monthly data for maximum Hs and companion Vw-in the regression analysis, to assess how adequate is the use of the return period value of the explanatory variable (in the regression analysis sense) to predict the return period value of the other variable and its uncertainty. This is to be carried out for each month of the year to estimate seasonal extreme (and companion) values associated with given return periods, so that they can be considered as a demand For the scope of this study, only models corresponding to the solid lines in Figures 1 and 2 are used (i.e., the whole set of available data at buoy 42,001 is considered), although projected metocean variables for other return periods are also used. The values for other return periods, used later in this study, for each month are listed in Appendix A in Tables A1 and A2 for Hs and V w , respectively. No details in this study are given about the parameters of the GEV models as were provided in previous studies (e.g., [28]), since attention is focused on proposing a simplified approach to correlate extreme Hs data and their companion V w data associated with given T -yr values. This is developed in the following section.

Simplified Approach for Return Period Correlated Values
The proposal to estimate the simultaneous occurrence of extreme metocean variables is based on a modified version of the classical linear regression technique. This consists in including estimated values from the time-dependent GEV models in the previous section-together with the simultaneously observed monthly data for maximum Hs and companion V w -in the regression analysis, to assess how adequate is the use of the return period value of the explanatory variable (in the regression analysis sense) to predict the return period value of the other variable and its uncertainty. This is to be carried out for each month of the year to estimate seasonal extreme (and companion) values associated with given return periods, so that they can be considered as a demand which is applied simultaneously to a maritime system (or other system). From this standpoint, significant wave heights and wind velocities can be considered as concurrent extreme metocean hazards accounting for seasonality, since they can be translated into load effects affecting a given maritime structure at the same time, at any desired season of the year. Moreover, such demands are usually related to given return periods for design purposes; consequently, it is convenient to establish them in terms of return periods. This is also the case for using contour maps in ocean engineering [20].
Before proceeding to introduce the proposal, the classical regression analysis using Hs as the explanatory variable and V w as the variable to be predicted is described as follows by using the observed data only. Buoy 42,001 is used to illustrate the procedure, but this can be extended to any buoy with simultaneous recorded metocean data. Hs is selected as the explanatory variable because it was used in a previous study [28] to investigate maxima. As noted before, V w is not necessarily the maximum, but the wind velocity observed simultaneously with the monthly maximum Hs. Naturally, it could be considered otherwise, e.g., maximum V w and associated simultaneously recorded companion Hs, but also other pairs of concurrent values of Hs and V w could be selected. This aspect is discussed later. From a designer perspective, it can be thought of as selecting an Hs value and using it to predict a V w value that is expected to occur simultaneously; then, both can be used to estimate the imposed (adding) effects of both metocean variables acting over a system, which capacity is to be designed to withstand such demand. Since both metocean variables are random, this can be expressed in mathematical terms [32,33] as: where the expected value of V w , given that Hs = hs, is a linear function of hs; b and α are constants to be determined by linear regression analysis. However, guideline expressions that relate wind velocities and wave heights, as well as fitted functions found in the literature [26], commonly have the following mathematical form: which can be linearized by taking the natural logarithm, leading to: where β = ln b; the notation in terms of the expected value is skipped for simplicity. Equations (7) and (9) are used in the linear regression analysis. For using Equation (9), the wind velocities and significant wave heights are first transformed into logarithmic space [33]. Figure 3a,b shows the linear regression analysis by using the nonlinear model (i.e., Equation (9)) and the linear model (i.e., Equation (7)), respectively, for October (around hurricane season); the residuals versus the fitted values of ln V w ( Figure 3c) and V w (Figure 3d) are also shown.
As mentioned before, the regression is first performed by using only the historical observed data, shown by circles in Figure 3a,b. It can be observed that the regression line (solid line) in both cases (i.e., Equations (9) and (7)) predicts increasing Vw for increasing Hs. This is consistent with what is reported in the literature. The 95% confidence intervals, which are measures of uncertainty for mean and future values of Vw [33,34], are also shown in Figure 3a,b (inner and outer dashed lines for mean  (9)) and (b) linear (Equation (7)) models for October. Residuals using (c) Equation (9) and (d) Equation (7).
Appl. Sci. 2020, 10, 4794 8 of 25 As mentioned before, the regression is first performed by using only the historical observed data, shown by circles in Figure 3a,b. It can be observed that the regression line (solid line) in both cases (i.e., Equations (9) and (7)) predicts increasing V w for increasing Hs. This is consistent with what is reported in the literature. The 95% confidence intervals, which are measures of uncertainty for mean and future values of V w [33,34], are also shown in Figure 3a,b (inner and outer dashed lines for mean and future values, respectively). Analysis of residuals versus predicted values indicates that the nonlinear model leads to a better scattering of the data along the horizontal axis (i.e., the model is uncorrelated) and better symmetry in the vertical axis (i.e., the zero mean assumption), indicating that the nonlinear model is better for the values observed in October. Both models roughly exhibit a normal distribution behavior in the residuals, as can be observed in Figure 4a for the model in logarithmic space and 4b for the linear model, where the normal probability papers are plotted using the cumulative probability given by i/(n + 1), where n is the number of samples (i.e., the number of residuals) and i corresponds to the number of a datapoint (a residual) in ascending order.    (9)) and (b) linear (Equation (7)) models for February. Residuals using (c) Equation (9) and (d) Equation (7).   (9) and (7), respectively. It is found in this study that correlation of significant weight height and wind velocity can be better represented (in terms of classic linear regression) by linear (Equation (7)) or power (Equation (9)) expressions, depending on the month of the year (e.g., hurricane season, cold front season or others). This is not clearly stated in other studies.
Nevertheless, it can be observed that both the linear or the nonlinear regression analysis lead to a fairly normal behavior of the residuals, and to a relatively uniform distribution of the variance around the mean, which could be associated with a homoscedastic behavior [33], although this varies depending on the month and type of regression considered (i.e., on using Equation (7) or Equation (9)). The error variance, σ e 2 is an indication of the uncertainty in the regression, and its squared root is known as the root mean squared error (RMSE). Overall, the residual analyses for all months indicate that the usual assumptions in the regression approach are met [34], and that Equation (7) and Equation (9) can be suitable alternatives. Although Equation (8) is frequently encountered in codified design and the literature [26] as indicated above, functional forms like those of Equation (7) are also found, as in the standard of the North Atlantic Treaty Organization, Standardized Wave and Wind Environments and Shipboard Reporting of Sea Conditions (NATO-STANAG 4194, see [26]), which can be expressed as: where the nomenclature used in the present study applies (for the NATO-STANAG 4194, b = 0 and α = 3.2). Equations with the functional form of Equations (8) and (10) are used for open sea and coastal waters; those employed for coastal waters are also linked to the fetch as an important variable [26]. In the present study, the fetch is not considered since the buoy used is located in open sea.   (9)) and (b) linear (Equation (7)) models for February. Residuals using (c) Equation (9) and (d) Equation (7).
Nevertheless, it can be observed that both the linear or the nonlinear regression analysis lead to a fairly normal behavior of the residuals, and to a relatively uniform distribution of the variance around the mean, which could be associated with a homoscedastic behavior [33], although this varies depending on the month and type of regression considered (i.e., on using Equation (7) or Equation (9)). The error variance, σe 2 is an indication of the uncertainty in the regression, and its squared root is known as the root mean squared error (RMSE). Overall, the residual analyses for all months indicate that the usual assumptions in the regression approach are met [34], and that Equation (7) and Equation (9) can be suitable alternatives. Although Equation (8) is frequently encountered in codified design and the literature [26] as indicated above, functional forms like those of Equation (7) are also found, as in the standard of the North Atlantic Treaty Organization, Standardized Wave and Wind Environments and Shipboard Reporting of Sea Conditions (NATO-STANAG 4194, see [26]), which can be expressed as:  (9)) and (b) linear (Equation (7)) models for February. Residuals using (c) Equation (9) and (d) Equation (7).
As shown in Figures 3 and 5, the data are approximately contained within the confidence intervals for future values in all cases. It is noticed that the confidence intervals for future values can be simplified by using Equation (11a,b), for the lower and upper interval, respectively.
where σ e = σ 2 e is the so-called RMSE which-advantageously-is part of the information that can be obtained when a regression analysis is performed; δ is a constant value to be discussed. If Equation (8) is to be considered to derive simplified expressions for the confidence intervals, expressions like Equation (11) can be written by including ±δ times the RMSE in Equation (9), except that the RMSE corresponds to the regression in logarithmic space. If the equation in the original space is to be used to compute the expected value (i.e., Equation (8)), it can be shown with some algebraic manipulations that the 95% confidence intervals for future values can be approximated by using Equation (12a,b): The value of δ is adequate to represent the 95% confidence intervals (denoted simply as confidence intervals from now on) for future values by adopting constant values approximately between 1.9 and 2.3, depending on the season considered and for both sets of equations (i.e., Equations (11) or (12)).
Exact confidence intervals (dashed lines) and approximated confidence intervals (Equations (11) or (12); dotted lines) are shown in Figure 6a-d for the four cases depicted in Figure 3a,b and Figure 5a,b, respectively; the original instead of the logarithmic scale is used for the cases corresponding to Equations (8) and (12). A value of δ equal to 2.3 is used for October (Figure 6a,b), while δ = 2.2 is employed for February (Figure 6c,d). It can be observed that Equations (11) and (12) are reasonable alternatives to represent the uncertainty in a simplified (practical) way. Even though simple, as far as the authors know, Equations (11) and (12) are not reported in the literature for extreme values accounting for seasonality; therefore, they are an important contribution of this study that include the uncertainty in correlated metocean variables for design purposes.  (11) and (12)) 95% confidence intervals, considering October (upper row; (a) and (b) for non-linear and linear cases, respectively) and February (lower row; (c) and (d) for non-linear and linear cases, respectively).
To investigate whether the previous proposal is consistent with both time-dependent GEV models in the previous section, estimated values of Hs and Vw associated with T-20, T-30, T-50, T-75 and T-100 are included in the scattergrams and regression analyses as can be observed in Figure 7, where the T-yr values for October are indicated by black filled circles with embedded perpendicular lines (lines to be explained below). Selection of these return periods is somewhat arbitrary but based on the usage of such return periods for design. Figure 7a,b correspond to the logarithmic and linear cases as before, with the regression line depicted with a solid line and confidence intervals with dashed lines (it is emphasized once more that they are derived by including the set of T-yr values from the GEV models in the regression). In fact, the regression lines and confidence intervals previously obtained (i.e., without including the black filled circles in the regression) are also plotted in Figure 7a,b using dashed-dotted and dotted lines, respectively. It can be observed that the regression lines are almost coincidental, that the estimated T-yr values are very near or over the regression line (especially for Figure 7b), and that the inclusion of T-yr values in the regression slightly decreases the uncertainty in terms of confidence intervals (i.e., the previously computed intervals are wider), which applies to both, Figure 7a,b.
Additionally, the amplitude between the lower and upper limits of the 95% confidence intervals from the GEV models for Hs and Vw are depicted with the perpendicular embedded lines in the T-yr values (i.e., embedded lines in the black filled circles) in Figure 7; vertical and horizontal embedded  (11) and (12)) 95% confidence intervals, considering October (upper row; (a) and (b) for non-linear and linear cases, respectively) and February (lower row; (c) and (d) for non-linear and linear cases, respectively).
The proposed approach to estimate concurrent metocean variables (Hs and V w in this case), is based on using the time-dependent GEV model for Hs obtained in the previous section, to compute a Hs value associated with a selected return period in years, T -yr (the subscript is used to denote the selected return period, e.g., T -30 corresponds to a 30-yr return period); the Hs for a given T -yr can then be used as input for Equations (8) and (10)- (12), to determine the corresponding expected companion V w for the same T -yr , including an uncertainty measure (i.e., the confidence intervals).
To investigate whether the previous proposal is consistent with both time-dependent GEV models in the previous section, estimated values of Hs and V w associated with T -20 , T -30 , T -50 , T -75 and T -100 are included in the scattergrams and regression analyses as can be observed in Figure 7, where the T -yr values for October are indicated by black filled circles with embedded perpendicular lines (lines to be explained below). Selection of these return periods is somewhat arbitrary but based on the usage of such return periods for design. A4l. Sci. 2020, 10, x FOR P9R REVIEW 12 of 25     Figure 6a,b and they are comparable, as can be observed in these figures, with slightly less wider confidence intervals in Figure 9; as in Figure 6 for October, δ = 2.3 is also used in Figure 9.  Figure 9b; all these equations are retrieved from [26], where details and proper citation can be found. Some of these equations are given as a function of Vw instead of Hs [26]; in such cases,  Figure 7a,b correspond to the logarithmic and linear cases as before, with the regression line depicted with a solid line and confidence intervals with dashed lines (it is emphasized once more that they are derived by including the set of T -yr values from the GEV models in the regression). In fact, the regression lines and confidence intervals previously obtained (i.e., without including the black filled circles in the regression) are also plotted in Figure 7a,b using dashed-dotted and dotted lines, respectively. It can be observed that the regression lines are almost coincidental, that the estimated T -yr values are very near or over the regression line (especially for Figure 7b), and that the inclusion of T -yr values in the regression slightly decreases the uncertainty in terms of confidence intervals (i.e., the previously computed intervals are wider), which applies to both, Figure 7a,b.
Additionally, the amplitude between the lower and upper limits of the 95% confidence intervals from the GEV models for Hs and V w are depicted with the perpendicular embedded lines in the T -yr values (i.e., embedded lines in the black filled circles) in Figure 7; vertical and horizontal embedded lines correspond to the confidence intervals for V w and Hs, respectively. The confidence intervals from the GEV models (the embedded lines in Figure 7) are obtained by using return period values (Equation (6)) and their variance, linked to the variance-covariance matrix and obtained through the delta method [29]. It can be observed in Figure 7a,b that confidence intervals from the GEV models surpass the limits of the confidence intervals for mean values (inner dashed and dotted lines) but are contained within the boundaries of those corresponding to future values (outer dashed and dotted lines) from the regression analysis; it can also be observed that intervals for Hs are wider than for V w . If it were of interest to delimit in a simplified way the confidence intervals from the GEV models for a given T -yr , Equations (11) and (12) could be an option by providing a suitable δ value. Neither this comparison between the uncertainty of the GEV models and the uncertainty from the linear regression, nor the idea of using Equations (11) and (12) to relate both uncertainties, are found in the literature. Therefore, this can be also considered as a contribution of the present study. Figure 8 shows the residuals and normal probability papers by including the T -yr values in the regression, where the residuals of T -yr values lie very close to or on the zero-mean line, and the probability papers preserve roughly similar trends than without considering the T -yr values and follow an approximately normal distribution.     Figure 6a,b and they are comparable, as can be observed in these figures, with slightly less wider confidence intervals in Figure 9; as in Figure 6 for October, δ = 2.3 is also used in Figure 9.  Figure 9b; all these equations are retrieved from [26], where details and proper citation can be found. Some of these equations are given as a function of Vw instead of Hs [26]; in such cases, equations are solved for Vw as a function of Hs. The PM, ITTC and West Shetland expressions have   Figure 6a,b and they are comparable, as can be observed in these figures, with slightly less wider confidence intervals in Figure 9; as in Figure 6 for October, δ = 2.3 is also used in Figure 9.  Figure 9b; all these equations are retrieved from [26], where details and proper citation can be found. Some of these equations are given as a function of V w instead of Hs [26]; in such cases, equations are solved for V w as a function of Hs. The PM, ITTC and West Shetland expressions have the functional form of Equation (8) with values of b equal to 6.376, 4.636 and 5.31, respectively, and values of α equal to 0.5, 0.709 and 0.603, respectively. Likewise, for the NATO-STANAG and 2013 Interim Guidelines, the corresponding functional form is given by Equation (10) with parameters b equal to 0 and 6.9 and α equal to 3.2 and 2.2, respectively. For the previous expressions Hs is given in m and V w in m/s. Note that the NATO-STANAG and 2013 Interim Guidelines equations are developed for coastal waters but included in Figure 9b for comparison. Figure 9 shows that the ITTC and the 2013 Interim Guidelines equations are similar in shape to those developed in this study, but under-or overestimating the predicted values for October the former and the latter, respectively. The other equations lead to not so dissimilar results as those reported here, and they are roughly in between the confidence intervals for future values.
To inspect how the proposed approach varies when considering seasonality, and how results compare with respect to the equations reported in [26] for other seasons, results like those shown in Figure 7 are also shown in Figures 10 and 11, but for the months of February and August, when the highest predicted T -yr values are obtained at cold front and hurricane seasons, respectively; in addition, results like those in Figure 9 are plotted in Figures 12 and 13. Unlike Figure 9, δ = 2.0 is used in Figures 12 and 13, since it leads to a better approximation of the exact confidence intervals. As mentioned before, δ can be varied to better represent the exact confidence intervals; nevertheless, the range of values is not large and δ can be roughly set between 1.9 and 2.3. equal to 0 and 6.9 and α equal to 3.2 and 2.2, respectively. For the previous expressions Hs is given in m and Vw in m/s. Note that the NATO-STANAG and 2013 Interim Guidelines equations are developed for coastal waters but included in Figure 9b for comparison. Figure 9 shows that the ITTC and the 2013 Interim Guidelines equations are similar in shape to those developed in this study, but underor overestimating the predicted values for October the former and the latter, respectively. The other equations lead to not so dissimilar results as those reported here, and they are roughly in between the confidence intervals for future values. To inspect how the proposed approach varies when considering seasonality, and how results compare with respect to the equations reported in [26] for other seasons, results like those shown in Figure 7 are also shown in Figures 10 and 11, but for the months of February and August, when the highest predicted T-yr values are obtained at cold front and hurricane seasons, respectively; in addition, results like those in Figure 9 are plotted in Figures 12 and 13. Unlike Figure 9, δ = 2.0 is used in Figures 12 and 13, since it leads to a better approximation of the exact confidence intervals. As mentioned before, δ can be varied to better represent the exact confidence intervals; nevertheless, the range of values is not large and δ can be roughly set between 1.9 and 2.3.  It can be observed from Figures 7, 10 and 11 that expected V w values as function of Hs depend on the season and month considered, that the uncertainty also varies for each case and that the regression is relatively irrespective of the selected scheme (linear or logarithmic space) for October and February, but it does vary when August is considered. Figure 11 for August corresponds to the largest estimated T -yr values for Hs (hurricane season) and, as shown in Figure 11 [26], at least for the case of V w and Hs and buoy 42001, but it is possibly the case for other metocean variables and buoys, in the Gulf of Mexico or anywhere else. Therefore, suitable sets of region-specific and site-specific equations should be developed to take into account seasonal effects.
Nonetheless, it is interesting to note that the equation developed in [26] for West Shetland is practically coincident with the regression line in Figure 13a at hurricane season. [28]. It can be observed from Figures 7, 10 and 11 that expected Vw values as function of Hs depend on the season and month considered, that the uncertainty also varies for each case and that the regression is relatively irrespective of the selected scheme (linear or logarithmic space) for October and February, but it does vary when August is considered. Figure 11 for August corresponds to the largest estimated T-yr values for Hs (hurricane season) and, as shown in Figure 11, the regression line predicts better estimates when the regression is performed in the logarithmic space; in fact, when the T-yr values are included in the regression (solid line), there is not significant deviation from the regression with the observed data only (dashed-dotted line); conversely, an important change in the slope of the regression line is observed for Figure 11b, implying better estimates of correlated values for given T-yr if the regression is carried out in the logarithmic space. Figures 9, 12 and 13 show that a single equation cannot capture the seasonality of correlated metocean variables, because it can overestimate or underestimate the correlated values when comparing against the equations compiled in [26], at least for the case of Vw and Hs and buoy 42001, For values associated with T -yr levels in Figure 11a (or in other figures), the uncertainty from the GEV models for each T -yr (i.e., the embedded perpendicular lines in the black filled circles) are a sort of envelope roughly covered by the confidence intervals for mean values in Figure 11a. Furthermore, such an envelope can be covered by using Equation (12) and a smaller value of δ (e.g., δ = 1) than that used for the confidence intervals of future values. Therefore, Equations (8) and (10) can be considered as possible design aids-simple expressions for design-like those contained in guidelines [26], but also accounting for the uncertainty in a simplified way (Equations (11) and (12)), either based on simplified confidence intervals for future values or confidence intervals within which the GEV uncertainty is confined. Since sets of such equations can be obtained for any desired season, the proposal here represents a simpler alternative than the more complex models reported in the introduction; consequently, they could be more attractive for designers since they are consistent with expressions recommended in guidelines.

Figure 12. Exact and approximated 95% confidence intervals for February considering the T-yr values: regression (a) in logarithmic space and (b) in in linear space. Legends indicate expressions referred in
Results for other months (not shown before) are reported in Figures A1-A9 in Appendix B. By observing Figures 7, 10 and 11 and Figures A1-A9, it can be concluded that the regression in the logarithmic or linear space could be both alternatives, with better results depending on the month considered. An exception is shown by Figure A8 for November, where neither a good fit is obtained for the logarithmic case ( Figure A8a) nor for the linear one ( Figure A8b) (also see Table 1). This is attributed to atypical simultaneous occurrence of larger Hs and smaller V w values, as can be observed in Figure A8, but also to a large dispersion of the data along the whole range of Hs values, leading to a regression line somewhat horizontal (and roughly around 12 m/s), implying no correlation of Hs and V w for November. It would be interesting to investigate in the future why this behavior is found specifically for November at buoy 42,001 in the Gulf of Mexico. For all the residual analyses shown in Figures 7, 10 and 11, and Figures A1-A9 for the observed data, mean values of these residuals are listed in Table A3 of Appendix C. In general, it can be observed in Table A3 that the zero mean assumption is met. Table 1. Fitted parameters of regression analysis for buoy 42,001 in this study considering seasonality.

Month For Equations (8) and (12)
For Equations (10) and (11) where β is the intercept obtained from the regression in the logarithmic space.
Results of the regression are listed for each month in Table 1 for buoy 42,001. Table 1 includes the regression parameters using only the observed data and the two approaches described above, i.e., the direct linear regression and the regression in the logarithmic space. It also includes σ e and relative errors of the predicted correlated V w in relation to those obtained from the GEV model for a selected return period (50 years) using: where the relative error is given in percentage, V bench is the computed T -50 value of V w using the GEV model and V est is obtained as proposed in this study, i.e., by using Hs for T -50 obtained from the GEV model as input value of Equations (8) and (10) (with regression parameters estimated from observed data only). Selection of T -50 is somewhat arbitrary and used to compare deviations of V w estimated with the GEV model with respect to this study's proposal, but any other T -yr could be considered. For a more detailed inspection of these deviations, Figures 7, 10 and 11 and those in the Appendix B can be used. Errors in Table 1 are an additional aid to compare which regression scheme leads to better results for any given month. The proposal in this study can be summarized in the following steps: (1) Select an Hs value associated with a T -yr and month of interest (e.g., see Table A1 in Appendix A) as input for Equations (8) or (10). (2) Select the season (month) for which the correlated values are desired, and the corresponding parameters (Table 1) to be also plugged into Equation (8) or (10), and compute V w . (3) Use Hs in step (1) and V w in step (2) as simultaneous demands associated with the selected T -yr and month acting over a maritime system to be designed. (4) If uncertainty is to be used by the designer, the corresponding σ e together with Equations (11) or (12) can be used to establish upper (or lower) limits by using a δ value (e.g., δ = 2 could be used to approximate the 95% confidence intervals for future values).

Discussion
For the proposal briefly outlined in the four steps of the previous section, it is noted that for step (1) the T -yr values of Hs are related to the GEV model, but any projected return level of the significant wave height available to the designer could be useful, at least as an approximation. Also, since for any given month selected, two options are available in step (2), namely, Equation (8) or Equation (10), the decision for the more suitable option can be based on the ε V reported in Table 1, as well as on Figures 7, 10 and 11 and Figures A1-A9. Justification for considering Hs and V w for T -yr in step (3) as simultaneously occurring metocean variables is based not only on the correlation of recorded (observed) data, but also in checking that the correlation of projected Hs and V w for T -yr , obtained from the time-dependent GEV models, are consistent with the trend of the observed data. Finally, the uncertainty obtained as per Equations (11) or (12) in step (4) approximately delimits all data points if δ = 2 is used (confidence intervals for future values); if the uncertainty from the GEV models is to be employed as reference (perhaps considering the perpendicular embedded lines as an envelope), lower δ values could be used, for instance δ = 1.5, which would be a conservative value to delimit the uncertainty given by the GEV models. Use of Equations (11) or (12) using δ values could be seen as factored demands like the ones used in codified design.
The proposed method allows the possibility of using other metocean variables, data from other buoys and hindcast data. It could use other probabilistic models. More than two variables could be considered if multilinear regression analysis is implemented. If additional recorded data becomes available for any month, the matrix approach in [35] can be implemented to update the parameters in Table 1.
Since the maximum Hs and companion V w may not yield the maximum demand in a given system-like the points in an environmental contour associated with a T -yr may not lead to the maximum demand in the system [20]-other pairs of simultaneously recorded data could be considered; one obvious combination is using maximum monthly V w and companion Hs, but also other data pairs could be used, for instance the second largest monthly Hs and companion V w , the second largest V w and companion Hs, the third largest monthly Hs and companion V w and so on. Each of the previous combinations would lead to pairs of T -yr values to which a system could be subjected to for estimating load effects, and the maximum demand could be used for design (in fact, different pairs of T -yr values could be more critical to different elements of the system). It could be investigated whether these pairs of T -yr values correspond to points of environmental contours proposed by other authors. Nonetheless, the idea is consistent with the approach of using multiple points from the contours to evaluate a system [20].

Conclusions
In this study, a simplified methodology to compute maximum significant wave heights and companion wind velocities associated with given return periods, accounting for seasonality, is presented. Simultaneous data from a buoy in the Gulf of Mexico are used. The approach is developed from projected return levels of significant wave heights based on a time-dependent GEV model and classical regression of the two metocean variables. A time-dependent GEV model for the companion wind velocities is also developed, to assess the adequacy of the method to predict the wind velocity as a function of significant wave height for a given return period.
It is found that the time window selected to estimate return period values of metocean variables can have an important impact in the predicted return levels for some seasons.
It is concluded that correlation of significant weight height and companion wind velocity can be adequately represented by linear or power equations, which could be easily implemented for design purposes, with different parameters to account for seasonality, but with the same functional form. Results are not very dissimilar with predictions from simplified equations in the literature or guidelines; however, it is found that a single equation with given parameters cannot capture the seasonality effects.
It is also found that the uncertainty in the predicted companion wind velocities as a function of significant wave height can be determined in a simplified way by using the root mean squared error from the regression analysis, expressed as a set of proposed equations to determine approximately (but closely) the 95% confidence intervals of future values. These same equations can be used to represent the envelope of the uncertainty estimated from the GEV models for different return periods.
It is considered that the proposed approach is a simple but adequate method to determine concurrent metocean hazards associated with given return periods, which could be imposed on a system to estimate the demand for design purposes, while also providing measures of uncertainty. The proposed expressions do not differ substantially to those provided in guidelines; thus, they could be amenable to designers, while also incorporating the seasonality effects in a simpler way as compared to other methods available in the literature. To the authors' knowledge, some of the findings in the present study-like the simplified proposal to include the uncertainty in the correlated metocean variables-are not reported elsewhere. Acknowledgments: We thank two anonymous reviewers for their comments, suggestions and constructive criticism. We also thank guess Editor Guido Ventura and the editorial team of Applied Sciences for their help in the editorial process.

Conflicts of Interest:
The authors declare no conflict of interest.                  Appendix C Table A3. Mean values of residuals for each month and regression scheme.

Month
Log. Space Lin. Space Appendix C