## 3 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:  
\[ \begin{array}{l}G(x)=\exp\left\{-\left[1+\xi\left(\frac{x-\mu}{\psi}\right) \right]_{+}^{-\frac{1}{\xi}}\right\}\rightarrow\xi\
eq 0\\ \\ G(x)=\exp\left\{-\exp\left[-\left(\frac{x-\mu}{\psi}\right)\right]\right\} \rightarrow\xi=0\end{array} \  
\[F_{t}(z)=\left\{ \begin{array}{l}\left\{\exp\left\{-\left[1+\xi(t)\left(\frac{z- \mu(t)}{\psi(t)}\right)\right]_{+}^{\frac{-1}{\xi(t)}}\right\}\right.\right.\\ \left.\left.\exp\left\{-\exp\left[-\left(\frac{z-\mu(t)}{\psi(t)}\right)\right] \right\}\right.\right.\end{array} \right. \tag{3}\]  
where \([a]_{+}=\max\) [\(a\),0]. The GEV parameters account for nonstationarity with the following expression:  
\[\theta=\left\{ \begin{array}{l}\mu(t)=\beta_{0}+\beta_{1}\cos(2\pi t)+\beta_ {2}sin(2\pi t)+\beta_{3}\cos(4\pi t)+\beta_{4}sin(4\pi t)\\ \psi(t)=\alpha_{0}+\alpha_{1}\cos(2\pi t)+\alpha_{2}sin(2\pi t)+\alpha_{3} \cos(4\pi t)+\alpha_{4}sin(4\pi t)\\ \xi(t)=\gamma_{0}+\gamma_{1}\cos(2\pi t)+\gamma_{2}sin(2\pi t)+\gamma_{3} \cos(4\pi t)+\gamma_{4}sin(4\pi t)\end{array} \right. \tag{4}\]  
where mean values are denoted with \(\beta_{0}\), \(\alpha_{0}\) and \(\gamma_{0}\), harmonic amplitudes by \(\beta_{i}\), \(\alpha_{i}\) and \(\gamma_{i}\) (for \(i=1\), 2) and subharmonic amplitudes by \(\beta_{i}\), \(\alpha_{i}\) and \(\gamma_{i}\) (for \(i=3\), 4); \(t\) is given in years. The maximum likelihood estimators for the location \(\mu\), scale \(\psi\) and shape \(\xi\) parameters vary during the year; the regression coefficients and the probability distribution parameters are determined through the method of maximum likelihood (MML), in which the likelihood function for the GEV model is given by:  
\[l(\theta|t_{i}x_{i})=-\sum_{i=1}^{m}\left\{\log\psi(t_{i}+\left(1+\frac{1}{ \xi(t_{i})}\right)\log\left[1+\xi(t_{i})\left(\frac{x_{i}-\mu(t_{i})}{\psi(t_{ i})}\right)\right]_{+}+\left[1+\xi(t_{i})\left(\frac{x_{i}-\mu(t_{i})}{\psi(t_{ i})}\right)\right]_{+}^{\frac{-1}{\xi(t_{i})}}\right\} \tag{5}\]  
As previously noted, the location \(\mu(t)>0\), scale \(\psi(t)>0\) and shape \(\xi\) (t) parameters are expressed in terms of harmonic functions, whose amplitudes are regression parameters to be mathematically estimated. 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,\theta)=X_{q}(\mu(t),\psi(t),\xi(t))=\left\{\begin{array}{l}\mu(t)- \frac{\psi(t)}{\xi(t)}\left\{1-\left[-\log(1-q)\right]^{-\xi(t  
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 \(H\)s.  
Figure 1 shows the maximum \(H\)s recorded monthly (depicted with dots) and the 30 yr return period values, \(T_{\text{-30}}\), estimated with the GEV model (depicted with a solid line).