Extreme Value Analysis of Ocean Still Water Levels along the USA East Coast—Case Study (Key West, Florida)

Abstract

This paper provides an Extreme Value Analysis (EVA) of the hourly still water level (SWL) record at Key West, Florida dating back to 1913 to understand the statistical likelihood of the combination of high predicted tides and the more dynamic influences (predominantly hurricane induced storm surges) that can drive ocean water levels higher at the coast. The impact of hurricane ‘Wilma’ in 2005 dominates the records, producing a super-elevation of the SWL above Mean Sea Level (MSL) of 1155 mm with an estimated return period of around 147 years. This paper explores the duality of increasing risks of oceanic inundation due to extreme events and increasing mean sea level over time, whilst also providing a range of recommendations for ensuring improved predictive model fitting and robustness of EVA for SWLs measured at tide gauges. When integrated with future IPCC AR6 sea level projections, the return level plots from the EVA provide decision makers with key guidance for design, strategic planning and climate change adaptation purposes at Key West, Florida.

1. Introduction

The physical, social, environmental and economic threat associated with rising global mean sea levels (GMSL) is now well established through the range of Intergovernmental Panel on Climate Change (IPCC) assessment reports (e.g., Ref. [1]) and associated technical studies and reports conducted more specifically for the United States context through the National Climate Assessment process [2,3], and other independent researchers (e.g., Refs. [4,5,6,7]). In 2021, the National Oceanic and Atmospheric Administration (NOAA) completed a detailed report examining sea level rise scenarios for the United States [8] that feed into the upcoming Fifth National Climate Assessment Report (or NCAR5).

In addition to these global and national assessments, the Office of Economic Demographic Research provides detailed annual assessments of flooding and sea level rise to the Florida Legislature consolidating all the literature and Federal Agency advice on the States exposure to these risks [9].

In essence, the risk associated with tidal (or oceanic) inundation at a given point in time is the combination of unrelated physical phenomena including

• high-intensity dynamic component (comprising the limit of any hurricane driven storm surge above the coincident tide condition);
• the prevailing tide coupled with any other omnipresent anomalies driven by prevailing oceanographic effects; and
• coincident, long timescale rise in relative MSL (due principally to climate change drivers and vertical land motion influences).

The State of Florida is especially vulnerable to climate change and concomitant sea level rise. Elevations in the coastal zone and southern peninsula are low and already subject to flooding during intense rainfall events, exceptional tides, and storm surge events [10]. On top of this, evidence highlights the northern foreshore of the Gulf Coast and along the east coast of the USA south of the Chesapeake Bay region being more exposed to the range of factors exacerbating threats from sea level rise than other coastlines at present [6].

This case study is based on Key West, Florida which examines extreme ocean water level phenomena along the south eastern coast of the USA. This paper is focused on understanding the peak heights (or extremes) measured above MSL resulting from the coincidence of predictable tides and the high-frequency dynamic influences, which are predominantly associated with Atlantic hurricane storm surges. The analysis benefits from the extensive hourly SWL recordings available from the Key West tide gauge facility dating back to 0600 h on 19 January 1913 (refer to Figure 1).

Extreme Value Analysis (EVA) techniques permit the fitting of continuous probability distribution functions which are used to predict rare (or extreme) phenomena often well beyond the boundary of the measured data [11]. The mathematical and statistical basis for EVA is well established in the literature (e.g., Refs. [12,13]) with ‘Block Maxima’ (BM) and ‘Peaks-Over-Threshold’ (POT) approaches amongst the most popular having been applied extensively across the engineering design, strategic planning, financial analysis, environmental and risk sectors.

The POT approach is more commonly applied for the consideration of environmental or natural phenomena as all independent events above a notional threshold are included. This is considered an advantage over the alternative BM approach which uses only a single maxima per equidistant time segments, ignoring other key rare events occurring within the same time segment [12]. For this reason, the current study uses the POT approach.

Notwithstanding the advantages of the POT approach to examine high-frequency (hourly) extreme ocean water level records, the utility of the approach is reliant on the selection of an appropriate threshold value above which to consider the peaks of declustered excesses via the fitting of the Generalized Pareto Distribution (GPD). The sensitivity of the threshold selected to drive the analysis presents as a rather subjective choice guided by expert judgement. Teena et al. (2012) [14] notes the choice of threshold is an important practical problem, which is mainly based on a compromise between bias and variance considerations. On one hand, the threshold must be high enough for the excess over the threshold to converge to the GPD whilst ensuring the sample size remains large enough to ensure sufficient data points for satisfactory determination of the GPD parameters [15].

A feature of the analysis contained herein is the treatment of the threshold selection based on several diagnostic treatments detailed in Section 2 which build on the existing literature providing guidance on best-practice applications for the POT technique (e.g., Refs. [11,16,17,18]) for natural phenomena with relevancy for coastal SWLs. These approaches have been largely automated using the ‘extRemes’ extension package in R [19,20].

This paper examines the high-intensity dynamic component associated with hurricane-driven storm surges which have dominated extreme hourly SWL measurements at the Key West tide gauge. Importantly, the paper also examines the influence of sea level rise over the data record and the impact of IPCC future projections associated with climate change. The paper also provides an assessment of the tide gauge record to consider the somewhat open academic question regarding the likelihood of increased storm intensity and winds to drive increased storm surges into the future associated with climate change.

The paper is structured with an explanation of the analytical methodologies applied and the data used (Section 2), leading into a presentation of the results of the current analysis (Section 3), a discussion section (Section 4), and, finally, conclusions (Section 5).

2. Materials and Methods

Various data sources and methodologies were applied in the study and are described in detail in the following sections.

2.1. Data Sources Used in the Study

A range of data sources were used to facilitate different parts of the analysis. Annual average time series data for Key West from the public archives of the Permanent Service for Mean Sea Level (PSMSL) [21,22] (Station ID = 188) were used in the MSL analysis spanning the period from 1913 to 2022 (inclusive). This data were used to determine the MSL trend over the course of the record.

High-frequency (hourly) observations spanning the timeframe from 19 January 1913 (0600 h) to 31 December 2019 (2300 h) were sourced from the University of Hawaii, Sea Level Centre (UHSLC) [23,24] (Station ID = 242). Only the ‘Research Quality’ data service available via the legacy data portal was sourced to undertake the EVA component. The residual of the hourly data up to 31 December 2022 (2300 h) was sourced from the NOAA’s Center for Operational Oceanographic Products and Services (CO-OPS) [25] (Station ID = 8724580 Key West, FL). Only the ‘verified’ form of the data was used and directly synchronized with the datum used for the hourly UHSLC data via direct and overlapping comparison with a common year of data (2019).

Sea level projections for a range of Shared Socioeconomic Pathway (SSP) scenarios used in AR6 were sourced directly from the IPCC providing global averaged sea level projections at decadal intervals between 2020 and 2150 [26,27,28] (refer to

Table 1. IPCC AR6 Global Sea Level Projections from 2020 to 2150 (in millimeters).

Year	SSP1-1.9	SSP1-2.6	SSP2-4.5	SSP3-7.0	SSP5-8.5
2020	0	0	0	0	0
2030	43	43	44	45	48
2040	79	87	93	96	107
2050	127	140	155	166	182
2060	161	183	213	233	260
2070	211	241	281	317	354
2080	252	289	354	409	459
2090	299	337	427	512	583
2100	335	387	507	630	716
2110	379	448	579	709	806
2120	417	496	653	819	930
2130	453	542	726	929	1050
2140	487	588	799	1037	1164
2150	520	632	870	1141	1272

Data sourced from [26,27,28] normalized to start date of 2020. Only median values advised.

for details).

For the preparation of various mapping products, hurricane track data were sourced from the US National Hurricane Center [29].

2.2. Methodology

The analytical steps involved in estimating extreme sea levels at Key West, Florida are detailed in the following sections. The application of fitting the GPD for EVA is based on the data conforming to the underpinning principles of stationarity and independence [12,13,30,31,32,33]. In practice, this requires conditioning the input data to ensure both the removal of any trend present (to satisfy stationarity criteria) and then declustering of input data to ensure multiple readings are not attributable to the same storm event (to satisfy independence criteria). These are the key steps prior to consideration of the EVA. All analysis and graphical outputs were developed by the author from a customized scripting code within the framework of the R Project for Statistical Computing [19].

Step 1: Determination of the MSL trend. This is a critical initial step as the fitting of a GPD function to estimate extreme values is based on the statistical principles of stationarity. All continuously recording tide gauges (such as Key West) measure water levels which are a function of all physical processes causing the water level to change, including but not limited to the melting of snow and ice reserves and thermosteric sea level change (resulting from climate change), dynamic climate mode (and other storm-related) influences and vertical land motions (VLM) at the tide gauge [6,34]. The isolation of the MSL trend so that it can be removed from the hourly water level record is no trivial exercise, but it is extremely important as the precision (or lack thereof) of the trend estimation can influence EVA [17].

Extensive time series analysis testing determined that data adaptive spectral techniques such as one-dimensional Singular Spectrum Analysis (SSA) provide greater utility in isolating the trend component with improved temporal precision for long tide gauge records [35,36]. Over recent years, this analysis procedure has been optimized for the benefit of sea level research in numerous updated regional mean sea level studies (e.g., Refs. [6,34,37]). In summary, the MSL trend is estimated from applying one-dimensional SSA, using a maximum (half-length) embedding dimension and aggregating components from the SSA decomposition which reside below the frequency band of 0.02 cycles per year (or above a 50-year period). For a more expansive discussion on the underpinning technique and parameterization of the SSA analysis, the reader is referred to Watson (2021) [6].

To apply the approach to the annual MSL data for Key West, the SSA procedure requires a complete time series. Therefore, the missing annual MSL data point for 1953 was filled using an iterative SSA procedure [38,39] within the R extension package ‘TrendSLR’ [19,40], which fills gaps based on the dominant spectral properties of the continuous sections of the original time series.

The annual MSL time series data from the PSMSL for Key West was initially synchronized with the datum used for the hourly UHSLC data via direct comparison to annual mean of hourly measurements. In this way, the extracted MSL can be directly applied to the hourly UHSLC dataset used for the EVA component.

Step 2: Detrending of the hourly tide gauge measurements. The comparatively smooth MSL trend derived in Step 1 needs to be converted from an annual time series to one that matches the hourly time series spanning 19 January 1913 (0600 h) to 31 December 2022 (2300 h). This is a relatively straightforward procedure whereby a cubic smoothing spline model is fitted to the isolated MSL trend from Step 1 to permit the prediction of MSL at each hourly time step over the time span of the annual data (May 1913 to May 2022). The cubic smoothing spline is extended linearly to encompass the residual hourly time steps beyond the periphery of the annual time series (approximately 6 months at either end). From here, the predicted hourly MSL trend is simply subtracted from the hourly tide gauge measurements to produce an hourly detrended (or stationary) dataset. The resultant time series represents hourly measurements above/below MSL.

Step 3: Declustering the hourly input data. Along with stationarity, statistical independence is the other key requirements for fitting a GPD to data for estimating extremes. Declustering procedures (i.e., making use only of the single highest exceedance within a cluster) are routinely employed in applications of the POT approach to avoid the effects of dependence [41]. The clustering influence of hourly tide gauge measurements have been examined previously [11,17], concluding no discernible influence from clustering on return interval water levels when the time span between specified events was set at >24 h. The ‘decluster’ function in the ‘extRemes’ package [19,20] has been used to decluster the detrended dataset (Step 2) by setting 25 h as the minimum span between successive peaks above the notional threshold [12,20]. Hourly measurements of extremes separated by fewer than 25 non-extremes are therefore considered to belong to the same cluster (or event).

Step 4: Extreme Value Analysis (EVA). Having addressed the conditions of stationarity (via Step 2) and dependance (via Step 3), the hourly input data file is now ready for EVA. The next key step for the application of the POT method is the selection of an appropriate threshold level above which to apply the GPD for EVA. The choice of the extreme threshold, where the GPD model provides a suitable approximation to the excess distribution, is critical in applications [42]. The choice of the extreme threshold above which to fit the GPD model involves balancing aspects of bias and variance [43].

On the one hand, if the threshold selected is too low, the GPD fit will be biased by the lower bound excesses which do not conform to the asymptotic tail laws underpinning the GPD. On the other hand, the choice of a threshold set too high will result in only a limited portion of the sample population from which to train the scale and shape parameters to fit a GPD model, resulting in excessive variance. There is an extensive body of literature dedicated to appropriate threshold selection. Coles (2001) [12] notes a range of common graphical diagnostics for threshold choice including the mean residual life (or mean excess) plot, threshold stability plot, and a suite of the usual distribution fit diagnostics (e.g., probability plots, quantile plots, return level plots, empirical and fitted density comparison).

It is acknowledged that some of these diagnostic tests and graphical results for threshold selection are challenging to interpret and require a degree of experience and skill [12,13,42,43,44]. Of these approaches, the mean excess over threshold approach [45] has been widely adopted throughout the literature and applied in this study. Simply described, the optimal advised threshold from mean excess over threshold plot can be determined as the lowest threshold for which the mean excess is approximately a straight line (within uncertainty bounds) [14,20].

Another factor which can affect the utility of EVA is the estimation method used to select the unknown parameters of the GPD tailed model fit. There are a wide variety of selection techniques which occupy the literature for this type of application. Yilmaz et al. (2021) [46] provided a detailed and wide-ranging testing regime of estimation methods for extreme value distributions encompassing classical techniques and Bayesian estimators applied to real-world data. Simulation results show that Bayes methods demonstrate better performance than the classical methods for estimating the unknown parameters. For this reason, a Bayesian estimation method was chosen from the available options in the ‘extRemes’ software package [19,20] despite the considerable computational disadvantage over the alternative methods available (MLE, GMLE and L-Moments).

Return interval plots and other diagnostic model fittings tools (such as quantile (or Q-Q) plots and root mean square error (RMSE)) were considered to ensure robust model fitting using the combination of selected threshold and estimation method. A detailed assessment of parameter sensitivity forms part of the Discussion section of the paper (refer to Section 4.2).

The ‘fevd’ function in the ‘extRemes’ package [19,20] was used to estimate the GPD model fit and to extrapolate empirical and model outputs (including confidence levels) for return interval plots used in the analysis and sensitivity checks.

3. Results

Results of the analysis undertaken are summarized under the following subsections for ease of presentation.

3.1. Hourly Tide Gauge Record

Figure 2 provides a pictorial summary of the detrended hourly tide gauge measurements spanning the timeframe from 19 January 1913 (0600 h) to 31 December 2022 (2300 h). Over this timeframe, there are 943,282 hourly records available for analysis with only 20,528 missing values (2.1%).

From this record, there are some 1146 measured hourly water levels above MSL that exceed the Highest Astronomical Tide (HAT) at the Key West tide gauge (540 mm). Following declustering, this reduces to 189 independent extreme events recorded above HAT (Highest Astronomical Tide).

The bottom panel of Figure 2 highlights the five largest independent events with the peak (1155 mm) recorded on 24 October 2005 (0700 h) resulting from hurricane ‘Wilma’ driven storm surges.

Table 2. Summary of 10 Largest Hourly Measurements above MSL at Key West.

Rank	Measurements above MSL (mm) 1	Date (Time, hrs)	Hurricane Event 2
1	1155	24 October 2005 (0700)	Wilma
2	989	10 September 2017 (1700)	Irma
3	907	8 September 1965 (1400)	Betsy
4	902	18 October 1944 (2200)	Unnamed
5	890	28 September 2022 (0300)	Ian
6	764	5 October 1933 (0300)	Unnamed
7	748	25 September 1998 (1800)	Georges
8	742	21 September 2005 (0200)	Rita
9	741	21 September 1948 (1700)	Unnamed
10	714	9 December 1946 (0700)	-

¹ Summary based on detrended and declustered results. Refer to Section 2.2 (Step 3) for further details. ² Hurricane event details sourced from National Hurricane Centre [29].

provides a summary of the 10 highest recorded extreme events at Key West, whilst Figure 3 provides a pictorial representation of the Atlantic hurricane storm tracks associated with the 5 highest recorded extreme events across the historical record. A total of 9 of the 10 highest extreme sea level events at Key West have occurred in September or October coinciding between the Atlantic hurricane season (1 June to 30 November [29]) and Florida’s ‘King Tide’ season (which occurs September through November [47]). Perhaps coincidentally, four of five extreme water level events recorded at Key West occurred with hurricanes tracking through a near identical crossover point (25° N/83° W), approximately 130 km to the northwest of Key West in the Gulf of Mexico.

3.2. Extreme Value Analysis

A central aspect of the POT/EVA technique involves the selection of an appropriate threshold above which extremes conform to the laws of the GPD whilst not being too high that such a small sample of the population results in excessive variance concerning the fitted GPD model. The mean excess plot (refer to Figure 4) was used as a tool to aid in the threshold selection. From this plot, the optimum threshold of 560 mm is readily identified as the lowest threshold for which the mean excess is approximately a straight line (within uncertainty bounds) by the orange-colored box. The selection of this threshold level is also prudent from a physical perspective, sitting marginally above the HAT for the Key West tide gauge (540 mm). The selection of threshold levels at or below the HAT increasingly bias the GPD model fit to excesses that do not conform to the statistical validation underpinning the GPD application.

Using the threshold value of 560 mm and a Bayesian estimation approach for the selection of the unknown GDP parameters, the EVA of the hourly height above MSL at the Key West tide gauge can be summarized in the return level plot in Figure 5.

Table 3. Predicted Extreme Hourly SWL above MSL (Key West, Florida).

Return Period (Years)	Height above MSL (mm)
1	576
2	613
5	674
10	733
20	807
50	937
100	1067
200	1236
500	1550
1000	1888

Central model estimates derived from Figure 5. Refer to Figure 5 for 95% CI.

summarizes extreme hourly heights above MSL for a range of return periods using the optimum fitted GPD function (Figure 5). Of interest here is the concave shape of the return level plot indicative of a heavy-tailed distribution. This suggests more randomness associated with the more extreme values, rather than the common bounded upper-tailed distributions which are convex in shape, asymptotically approaching some naturally constrained upper bound. The sensitivity testing of the GPD model fit for threshold and parameter estimation are explored more fully in Discussion (refer to Section 4.2) but confirm the high confidence in the fitted model for EVA at Key West.

The extreme value analysis highlights the impact of storm surges driven by hurricanes, no more so than hurricane ‘Wilma’, which dwarfs all other events recorded at this site. This event produced a super-elevation of the SWL above MSL of 1155 mm on 24 October 2005, with an estimated return period of around 147 years. The second largest event on record occurred on the 10 September 2017 in response to the storm surge driven by hurricane ‘Irma’. The super-elevation of the SWL above MSL of 989 mm equates to an estimated return period of around 67 years. For comparison, the HAT at Key West measures ≈ 540 mm above MSL [25].

Importantly, the analysis undertaken shows a relative rise in MSL over the period from 1913 to 2023 of approximately 280 mm. The effect of this sea level rise over that period is graphically illustrated in Figure 6, highlighted by the fact that a 100-year ARI (Average Recurrence Interval) extreme storm surge above MSL (1067 mm) in 1913 would only require a 17-year ARI event in 2023 to attain the same static water level. The influence of future projections of sea level rise are further considered in Discussion (refer to Section 4.4).

4. Discussion

This study gives rise to various discussion points, highlighted in the following sections.

4.1. Is There Any Evidence to Suggest Extremes Have Been Increasing over Time?

It remains a somewhat open academic question regarding the likelihood of increased storm intensity and winds to drive increased storm surges into the future associated with climate change. Reed et al. (2022) [48] note that quantifying the impacts of the increasing sea surface temperatures on hurricanes and tropical cyclones globally remains a scientific challenge given the competing climate effects on other environmental parameters important for storm genesis and development, including wind shear and atmospheric stability [49].

Figure 7 provides a summary of the slope of a linear regression analysis of detrended and declustered extremes above a range of threshold levels for confidence levels including 80%, 90%, 95% and 99%. Considering the optimal threshold used for the EVA analysis (560 mm), there is a small positive slope of 0.2 mm/year based on 116 extreme events; however, for each of the confidence levels considered, the slope is not statistically different to zero. At a threshold level of 600 mm, the positive slope increases to 0.5 mm/year for 63 extreme events, but similarly, at all confidence levels considered, the slope is not statistically different from zero. At a threshold level of 700 mm, the positive slope increases to 1.8 mm/year for mere 10 extreme events on record. At the 80% confidence level, the slope is marginally statistically above zero (inferring a positive trend). Overall, there is no compelling statistical evidence of independent extreme events increasing over time at Key West (yet), though it is worth noting some 4 of the top 10 (and 7 of the top 20) extreme events have occurred after 2004, within the last 20 years.

Should this latter period prevalence of extreme events continue over the forthcoming decade, it is likely a small statistically significant positive linear trend will start to emerge.

4.2. Sensitivity Analysis Testing for EVA

Figure 8 summarizes the sensitivity analysis to select the optimum fitted GPD for extreme value prediction based on the four separate parameter optimization approaches (i.e., MLE, GMLE, L-Moments and Bayesian). Results are presented for threshold values ranging from 300 to 700 mm at 5 mm increments.

This detailed summary plot highlights the sensitivity of both threshold selection and parameter estimation method on the fitted GPD and resultant predictive capability for extreme event return periods. In essence, the POT/GPD method does require some necessary expert judgement and diagnostic assessment to improve the reliability and confidence in the model fit. Hawkes et al. (2008) [50] notes the importance of the choice of an appropriate distribution function which should be guided not only by a goodness-of-fit test but also by the robustness of the fit.

The following example highlights the importance of balancing these concepts to select an appropriate fitted model. Arns et al. (2013) provides a range of best-practice advice for EVA of ocean water level records based on extensive data assessment of tide gauge records around the German Bight and various international data sets. The strong recommendation was to use the POT/GPD approach with a selected threshold level of the 99.7th percentile of high water peaks, concluding this resulted in the most stable return water level estimates. This current study finds that the use of the 99.7th percentile of high water peaks to estimate the threshold level for EVA at Key West would produce near identical return water level estimates to those advised. However, the threshold issue alone does not of itself automatically result in the optimum fitted model in the current study.

Figure 9 summarizes the return level plots using the recommended threshold level of 560 mm, suggested by the mean excess plot (Figure 4). The results provide a clear visual perspective on the goodness-of-fit for the fitted GPD model under all parameter estimate methods within uncertainty bounds. In this example, the final choice of the Bayesian approach to estimate the GPD model parameters is confirmed in the bottom right panel of Figure 9 which highlights both the best model fit to the top 10 extreme values recorded and satisfactory error estimates. The improved result using the Bayesian estimation method accords with the findings of Yilmaz et al. (2021) [46].

4.3. The Shape of the EVA Return Level Plot

The concave shape of the return level plot (Figure 5) is indicative of a heavy-tailed distribution, suggesting more randomness associated with the more extreme values. This tends to contrast to the shape expected in natural system phenomena which one might anticipate would follow a bounded upper-tailed distribution which are convex in shape, and most perceived to be asymptotically constrained by natural forcings. For example, wave formations formed by wind are commonly constrained by the limit of measured wind speeds and available fetch distances across specific ocean basins.

However, in the case of the super-elevation of the water surface above MSL at Key West, there are a range of random factors that could contribute to raising the ocean water surface due to a high-intensity-hurricane-driven event. These include, but are not necessarily limited to, the associated storm surge from winds and atmospheric forcings, the path of the hurricane, pre-existing currents and thermosteric anomalies in the Gulf of Mexico, and any additive impacts from concomitant rainfall intensity, duration, and scale.

All the aforementioned factors are most likely to increase in intensity with climate change [46]. Although this study finds no statistical evidence (yet) of an increase in extreme water levels over time (refer to Section 4.1), it is likely that a small positive trend will emerge should the prevalence of extreme events measured over the past two decades continue over the forthcoming decade. Noting 7 of the top 20 extreme measured water level events have occurred over the past 17 years, it is not unexpected that the extreme predictions are well modelled by a heavy tailed distribution. Unpacking these issues further is beyond the scope of the current study, but it would be a very interesting and valuable academic endeavor.

4.4. Influence of Future Sea Level Projections on Predictions of Extreme Water Levels

The impact of rising mean sea levels projected over the course of the next century (and beyond) will be profound in exacerbating the already devastating consequences associated with the super-elevation of the ocean water level surface due to extreme phenomena. For example, an extreme hourly height above MSL with an annual expected recurrence (576 mm, refer to Table 3) is anticipated to be eclipsed by global mean sea level alone by around 2140 under a SSP1-2.6 projection scenario and by as early as 2090 under a SSP5-8.5 projection scenario (refer to Table 1). Considered from another perspective, the height above MSL reached by a 500-year extreme water level event in 2020 (1550 mm) would only require an extreme event with an estimated return period a little over 20 years in 2100 under global mean sea level rise associated with a SSP5-8.5 projection scenario.

The influence of these projections could be foreshortened by the additional hazard presented by vertical land motion. For example, it is estimated at the Key West tidal facility that land subsidence is occurring at a rate of around 1.2 mm/year [6] which would further exacerbate the influence of global mean sea level rise in the local context.

4.5. Recommendations for EVA Applied to Ocean Water Level Records

There is an extensive body of literature dedicated to the statistical application and utility of extreme value analysis across the varied disciplines in which it has been applied (e.g., Refs. [11,12,13,14,15,16,17,18]). There are some key suggestions for improving the utility and robustness of EVA applied to still water levels from long tide gauge records which can be broadly summarized as follows:

• use hourly, quality controlled tide gauge data;
• detrend and decluster the hourly data inputs (to attain stationarity and independence requirements for EVA);
• apply the POT approach with the GPD model fit;
• threshold selection recommended above HAT;
• optimize threshold selection through the use of tools such as mean excess plot;
• test a range of thresholds to better understand key sensitivities (if necessary);
• use Bayesian method to optimize parameter estimation; and
• confirm optimum threshold and parameter estimation method through visual diagnostic means including return level plots and Q-Q plots to observe the performance of the model fit within uncertainty bounds to the more extreme values.

4.6. Limitations of the Analysis and Results

It should, however, be noted that extreme value analysis involves a relatively small proportion of the overall data set which must adhere to quite specific mathematical foundations (i.e., stationarity and independence). As a result, the fitting of continuous probability distribution functions can be very sensitive to a small number of very extreme events resulting in quite large confidence intervals for increasing event rarity as evidenced in the return level plot for Key West (Figure 5). The key results of the analysis provide return levels for extreme water levels above mean sea level that are the combination of the tidal influence at Key West coupled with the meteorological forcings inherent in the regional context. Therefore, the results are not directly transferable to another location.

It is also noted that a small statistically significant positive trend in extremes might begin to emerge within the next decade or so. Should this eventuate, the extreme value analysis and return level plots espoused in this paper would no longer be valid (with stationarity principles violated) and the analysis would have to be reconsidered.

5. Conclusions

The extensive data repository of measured hourly still water levels at the Key West tidal facility dating back to January 1913 provides an outstanding resource to consider both high-frequency extreme water levels and low-frequency sea level rise in Florida.

The largest hourly measurement above MSL (1155 mm) at Key West was recorded on 24 October 2005 (0700 h) resulting from hurricane ‘Wilma’ driven storm surges. The recurrence interval of this event is estimated at ≈147 years. ‘Wilma’ formed and became an extremely intense hurricane over the northwestern Caribbean Sea. This system produced the all-time lowest central pressure for an Atlantic basin hurricane (≈882 mb), devastating the northeastern Yucatan Peninsula and inflicting extensive damage over southern Florida. Within 24 h, the event intensified from a 60 kt tropical storm to a category 5 hurricane, an unprecedented event for an Atlantic tropical cyclone [51]. The super-elevation of the SWL above MSL for the hurricane ‘Wilma’ event was some 166 mm higher than the second largest extreme driven by hurricane ‘Irma’ which was recorded on 10 September 2017 (1700 h).

Importantly, over the 109-year station record, there has been a relative increase in MSL of approximately 280 mm (due to climate change increases and subsidence at the tide gauge). The exacerbation of risk from inundation due to rising MSL over time is highlighted by the fact that a 100-year ARI extreme storm surge above MSL (1067 mm) in 1913 would only require a 17-year ARI event in 2023 to attain the same static water level. IPCC projected rises in global MSL from the 2020–2150 range from around 520 mm to 1272 mm depending on the modelled scenario (refer to Table 1).

It is perhaps unsurprising that 9 of the 10 highest extreme sea level events at Key West occurred in September or October, coinciding with both the Atlantic hurricane season and Florida’s ‘King Tide’ season. Interestingly, four of the five highest surge-producing events have resulted from hurricane’s tracking through a common point, approximately 130 km to the northwest of the tide gauge in the Gulf of Mexico. An extension to the current paper might investigate this further along with the role that climate change might play with influencing (or shifting) track paths of extreme events into the future and how this might impact extreme water levels not only at Key West, but also regionally.

Based on the data analysis herein, there is no statistical evidence of independent extreme events increasing over time at Key West (yet), noting 7 of the top 20 extreme events have occurred since 2004. It is, however, suggested that should this prevalence of extreme events continue over the forthcoming decade, it is likely that a small statistically significant positive trend in extremes will begin to emerge.

This paper provides an updated, state-of-the-art guide to robust EVA applied to long hourly SWLs recorded at a tide gauge for application at any location (refer to Section 4.5). The resultant hourly sea level extremes above MSL for a range of return periods (Table 3), coupled with future sea level projections (Table 1), provide sound resources for coastal design, planning and sea level adaptation purposes at Key West, Florida.