Empirical and numerical data were used for the theoretical analysis presented here. Both sources provide data on the eastward and northward components of the ocean currents in the study area for different depths and different temporal resolutions. The Canek project, which has carried out similar measurements in the past (Chávez et al., 2003), was responsible for the measurements of the current in the Cozumel Channel. The Canek research project, also known as the Estudio de la circulación y elintercambio a través del Canal de Yucatán (Study of circulation and exchange through the Yucatan Canal) has been coordinated by the Centro de Investigación Científica y Educación Superior de Ensenada since its foundation in 1996. The data were obtained using a stationary, long-range acoustic Doppler current profiler at N20°32.218′ W087°02.738′ [see Figure 1B] anchored at a depth of approximately 400m, and measuring every half hour, from 9^th April 2009 to 14^th May 2011. The data on depth were recorded in 16 cells, with the shallowest cell at a depth of 49m. For the numerical data, the HYCOM was chosen because of its good temporal range and resolution, and excellent spatial resolution, compared to other products. In this study, the data of the reanalysis model HYCOM + NCODA GOMu0.04 experiment 50.1 are used, which are publicly available in https://www.hycom.org/data/gomu0pt04/expt-50pt1. The model provides the current components at 40 depths for the Mexican Caribbean (among other regions), covering 1^st January 1993 to 31^st December 2012, at a temporal resolution of three hours and a spatial resolution of 0.04 in both eastern and northern directions. Numerical current data are reported at a depth of 50m while the empirical data describe the current at 49m.

Interpolation of HYCOM data to match the Canek data was carried out using the griddata function, available in the SciPy-module (version 1.6.1) for Python 3 (Virtanen et al., 2020). The four nodes of the numerical model were used as input, which surround the location of the measured field data. Due to the different sampling frequencies of the data sources, the data with higher frequency (i.e., the empirical data) had to be reduced. The data provided by the Canek project were sampled every 30 min and every hour, depending on the date. As the HYCOM numerical data reports the instantaneous value every three-hours, the empirical data were reduced by discarding every time step which is not available in the numerical data.

A linear quantile regression was performed on the current speed, using the quantreg model, as provided by the statsmodels module (version 0.12.2) for Python 3 (Seabold and Perktold, 2010). The linear regression was estimated for the empirical data proportional to the numerical data with the intercept set to zero.

To estimate the error produced by the numerical data, a simplified EVA was performed for both data sets, the empirical data in its original form and the numerical data in its adjusted form but reduced to the temporal range of the empirical data set. The analysis is described in detail in section 2.3. However, due to the low number of observations available in both sets, the methodology had to be modified. As a threshold, the 0.5^th-quantile was used in contrast to the suggested automated threshold selection. However, the same range of possible thresholds was used to estimate the confidence interval. The (signed) relative error between empirical and numerical data is defined as

where u_e is the empirical current speed and u_m is the speed as predicted by the numerical model. Besides the mentioned modules for Python 3 (Van Rossum and Drake, 2009), substantial parts of the data processing have been carried out with the NumPy-module in version 1.20.1 (Harris et al., 2020) and the pandas-module in version 1.2.2 (Wes McKinney, 2010).

The methodology used assumes that for a random variable (x) the excess over a suitable threshold (u) can be modelled by a GPD. Liang et al. (2019) define the cumulative density function of the GPD as

Where x represents the random variable, u the threshold, ξ the shape parameter, and σ the scale parameter.

The procedure can be summarized as follows, where the automated threshold selection method is based on the work of Thompson et al. (2009):

To consider the phenomenon as random, the realization of each variable should be independent. However, with the temporal resolution provided, the data analysed in this study is not random. To select only values independent of temporally close values (later called peaks), a moving time window was used, as suggested in Solari et al. (2017). The time window is of fixed length, depending on the variable type, and moves consecutively through the time series. If the value in the centre of the time window is the maximum of that time window, this value is selected as an independent peak.

Outliers may be present in the set of selected peaks, which would alter the final excess model. For the automated and semi-automated detection of outliers, a great variety of methods are available (Hodge and Austin, 2004). One of the simplest methods suitable for univariate data is based on quartiles and presented in Laurikkala et al. (2000). The authors define an upper (u_u ) and lower threshold (u_l ), beyond which all peaks are considered as outliers and are consequently discarded. Both thresholds are defined by

where q₁ is the first quartile (25^th percentile) and q₃ is the third quartile (75^th percentile).

From the previously selected peaks, potential thresholds are selected, as suggested in Thompson et al. (2009). The potential thresholds are equally spaced between 25^th and 98^th percentile. If less than 100 peaks are found above the 98^th percentile, the 100^th largest peak is selected as the upper limit of the range for potential thresholds.

For each threshold, all peaks x_i > u_j are selected, and a GPD is fitted through those peaks. The shape (ξ_j ) and scale (σ_j ) parameters of the GPD are determined by the function genpareto.fit, which is part of the SciPy.stats-package. The location parameter is held fixed to the corresponding threshold u. The reparameterised scale parameter, which is defined by

should be constant above a suitable threshold, following Coles (2001). This relationship was extended by Thompson et al. (2009) by fitting a normal distribution with a mean of zero through the difference of the reparameterised scale parameter for the current and all greater thresholds. This difference is defined by

The first threshold for which the corresponding normal distribution has a p-value ≥ 0.05 is selected for calculating the return levels. As a test for normality, the Kolmogorov-Smirnov test is used, as implemented in the ks_1samp function of the SciPy.stats-package. The return level X_m (Coles, 2001) can be calculated as

The average number of peaks m during a return period (T_B ) is defined by

where n_p is the total number of peaks and n_y the number of years for which data is available.

The exceedance probability of threshold ζ ^ u , the complete Variance-Covariance Matrix V and the variance of return level Var   ( x ^ m ) are estimated, as stated in the following equations (Coles, 2001), where the values with a hat indicate the estimation of the corresponding value.

where n_pot is the number of peaks over threshold

with

The empirical and numerical data used have different temporal resolutions. Therefore, the empirical data from the Canek project were downsampled by discarding all time steps that are unavailable in the data provided by the numerical model. In Figure 2 the unadjusted numerical data show a clear bias towards overestimation. The numerical data were adjusted by the linear regression model, which was estimated with a quantile regression using the 0.5^th-quantile, and is defined by

where u^′_m is the adjusted current speed. The adjusted numerical data reflect the empirical data much better. The effect of the model adjustment is strongly reflected by the mean relative error, that is reduced from −0.255 to 0.006. The mean absolute relative error of 0.288 and the root mean squared relative error of 0.365 are reduced to 0.153 and 0.206, respectively. Of greater concern, however, is the missing tail of the probability distribution of both numerical data in Figurs 2A, as these are of great importance for EVAs. Especially in the adjusted data, this leads to a pronounced underestimation of higher current speeds (i.e., rare events) as seen in the deviation from the diagonal in Figure 2B.

To quantify the effect of the missing tail on the extreme value estimations, a simplified EVA was performed. Due to the short time range, the length of the time window was reduced to 7 days (i.e., 56 observations). Additionally, the 0.5^th-quantile was selected as the threshold rather than the proposed automated threshold selection. As expected from the results presented in Figure 2B, the adjusted numerical model shows an underestimation of extreme values, as can be seen in Figure 3A. Nevertheless, for rare events (return period > 10 years) the relative error converges to a value just below 0.22 (see Figure 3B). The large 95% confidence interval (CI_95%) in Figure 3A is the result of the short temporal coverage of the data used for the analysis. It should be noted that the CI in Figure 3B is not CI_95%; it is the maximum error estimated by the upper and lower bounds of the CI_95% in Figure 3A.

The methodology was applied to several nodes in the Mexican Caribbean, shown as grey dots in Figure 4. It should be noted that not every node contains data on the current, as some are on land, or in waters of less than 50m depth. The four nodes marked in red were selected as the results suggest that it is possible to obtain a different behaviour with respect to the GPD fit. The node at position P1 (20.520° N, 86.600° W) is where the current is most concentrated off the east coast of Cozumel. The node at position P2 (20.640° N, 86.960° W), is in the Cozumel Channel, near a possible site for the installation of ocean current turbines [see Alcérreca-Huerta et al. (2019)]. That at P3 (21.040° N, 86.560° W) is in the wake of the Cozumel Channel, off the coast of Cancun, and the node at P4 (21.800° N, 86.480° W) is in the Yucatan current northeast of Cancun.

To determine the optimal length of the time window, the number of peaks identified in the windows was analysed for the nodes at P1 to P4. In Figure 5A, there is a steady fall in the number of peaks, but it remains above the critical number of 200. At a length of 25 days, the number of peaks for all four nodes is just below 250. The relative difference between the number of peaks and the length of time window (see Figure 5B) shows a decreasing trend, as the length of the time window increases. From a 21 day length, the relative difference is less than 10%, dipping briefly below the 5% mark at a length of 23 days. To spread the number of peaks evenly within the time range analysed, and to avoid having too few peaks, a time window of 23 days in length was chosen.

Figure 6 shows the statistical data of the GPD fit for each node. While the north, east, and southern boundaries of the domain are determined by the node selection, the western boundary is a feature of the numerical data generated by HYCOM for this site. The number of identified peaks seems to be similar in the study region (see Figure 6A), with a slight decrease towards deeper waters.

Figures 6B, C) show the selected threshold for each node and the corresponding p-value for the automatic threshold selection, respectively. The value of selected thresholds tend to increase in the centre of the channel, and in the stream close to the east coast of Cozumel Island that extends northward, along the Cancun coast. This is expected, since the current becomes more intense at these locations. As it was possible to find a suitable threshold for all nodes with information on the current velocity, the p-value is over 0.05 in significance, although some inconsistencies of above 0.1, and even 0.15, are found throughout the domain.

There is no clear trend in the shape parameter of the fitted GPD in Figure 6D. However, it was estimated to be negative for all nodes, producing a bounded GPD. For a few nodes at the northwestern boundary, the shape parameter was estimated to be very close to zero. The scale parameter in Figure 6E indicates a slight increase off Cozumel Island and at the northeastern boundary, which leads to a thicker tail for the GPD in those regions (i.e., increased return levels).

Figure 6F shows the number of peaks above the threshold which lie within the estimated CI_95%. For nearly all the nodes, the estimated CI_95% covers 100% of the numerical observations. As is to be expected, not all the observations are within the CI_95% for all the nodes. However, the number of nodes for which some observations are outside the CI_95% is small, while the minimum share within the analysed region is still above 90%. Despite this apparent overestimation of the CI_95%, this suggests that the methodology of GPD fit together with the estimation of the CI_95% is suitable and the results of the GPD for the given input data is reliable.

Return levels for the selected return periods, on the corresponding lower and upper boundaries of the CI_95% are shown in Figure 7. The expected return level (central column) increases in the channel and the main current, which extends northwards from the east of the Cozumel Island. This trend is further pronounced in the case of the CI_95% upper boundary (right column of Figure 7), which is in agreement with the results in Figures 6B, E. The region with higher shape parameters at the northwest edge of the domain (see Figure 6D) is not noticeable in the estimated return level (centre column in Figure 7). However, in case of the CI95%-limits, that region stands out with lower return level for the lower bound of the CI_95% and higher return levels for the upper limit, suggesting a much higher uncertainty. The distribution over the rest of the domain is as expected, see Figure 6.

The parameters for the GPD excess model for the four nodes seen in Figure 4 are summarized in Table 1. Except for the node at position P4, the shape parameters are negative, with all of CI_95% below zero. Compared to the standard error, the shape parameter at P4 is small, giving a CI_95% closely centred around zero. However, this result could be due to an error in the numerical model, as mentioned above. The highest scale parameter is found at P1, which produces a thicker tail to the probability distribution. Nevertheless, the bound nature of the excess model (due to the negative shape parameter) prevents high return levels for this node. The number of peaks found for each node is similar, just above the critical threshold of 200. Slightly more than 100 peaks were found above the selected threshold. The number of peaks, and peaks above the threshold, suggests that the selected time span of 20 years is a bit short, but still sufficient to perform the EVA.

In Figure 8, the peaks, POT, and the thresholds are shown for the four nodes. None of these nodes have a cluster of peaks (or lack thereof), suggesting that a 23 day time window is sufficient. The distributions of peaks, together with the filtered outliers, are shown for each node in Figure 9. The nodes at P1 to P3 show a standard distribution of peaks. The node at P4 has a multi-modal distribution, suggesting an error, and that the conclusions drawn from the data might not be reliable. No outliers at the upper end were found for the node at P2, whereas at P1 and P4 there were one each, and at P3, two. At the lower end, a few outliers were also detected and filtered out, but due to the nature of POT methods, these tend to have no significant effect on the outcome.

Figures 10–13 present the corresponding diagnostic plots for the GPD excess model. For a detailed interpretation of this type of plot, the reader is referred to Coles (2001). Despite slight deviations in the diagnostic plot for the node at P1 in Figure 10, and especially the q-q plot in Figure 10B, 100% of the empirical POT still lie within the CI_95%, as seen in Figure 10C and tab. 1. Both plots suggest an overestimation of the GPD model. There are few peaks, especially visible in the density plot (Figure 10D). However, the bound excess model seems to give a good fit for the underlying numerical data. The diagnostic plots for the node at P2 (Figure 11) show some deviations between the numerical data and GPD excess model in the p-p plot (Figure 11A). Around the 0.6 mark, the GPD excess model shows a slight overestimation. This deviation is also visible in the q-q plot (Figure 11B) at speeds of about 1.45ms−1 and in the return level-plot (Figure 11C) at the same speeds. Despite these inconsistencies, the GPD excess model is a fits the data well.

The p-p plot presents some discrepancies at the 60% percentile at P3 (Figure 12A). The excess model also differs from the numerical data for higher speeds, as seen in Figures 12B, C. However, all the observations are within the estimated CI_95% of the return level, suggesting that the GPD excess model application is reliable.

The plots in Figure 13 bring into doubt whether this excess model can be used to reliably estimate the extreme values of the node at P4. Although the CI_95% includes all the numerical observations, the p-p plot (Figure 13A) and especially the q-q plot (Figure 13B) looks unusual. A slight s-shaped deviation is present, with considerable inconsistencies above 1.2ms⁻¹ in the q-q plot. Additionally, and as observed in Figure 7, the CI_95% in Figure 13C is quite large, while the density plot in Figure 13D shows a reasonable fit to the data.