Statistical estimation of spatial wave extremes for tropical cyclones from small data samples: validation of the STM-E approach using long-term synthetic cyclone data for the Caribbean Sea

Occurrences of tropical cyclones at a location are rare, and for many locations, only short periods of observations or hindcasts are available. Hence, estimation of return values (corresponding to a period considerably longer than that for which data is available) for cyclone-induced significant wave height (SWH) from small samples is challenging. The STM-E (space-time maximum and exposure) model was developed to provide reduced bias in estimates of return values compared to competitor approaches in such situations, and realistic estimates of return value uncertainty. STM-E exploits data from a spatial neighbourhood satisfying certain conditions, rather than data from a single location, for return value estimation. This article provides critical assessment of the STM-E model for tropical cyclones in the Caribbean Sea near Guadeloupe for which a large database of synthetic cyclones is available, corresponding to more than 3,000 years of observation. Results indicate that STM-E yields values for the 500-year return value of SWH and its variability, estimated from 200 years of cyclone data, consistent with direct empirical estimates obtained by sampling 500 years of data from the full synthetic cyclone database; similar results were found for estimation of the 100-year return value from samples corresponding to approximately 50 years of data. In general, STM-E also provides reduced bias and more realistic uncertainty estimates for return values relative to single location analysis.

tions regarding the physical processes involved. Statistical estimation is problematic, since inferences must be made concerning extreme quantiles of the distribution of quantities such as SWH, using a limited set of data.
Objective and layout 55 In the present work, we aim to tackle the problem of realistic return value estimation for small samples of tropical cyclones using a recently-developed procedure named STM-E, which has already been successfully applied in regions exposed to tropical cyclones near Japan (Wada et al. 2018) and in Gulf of Mexico (Wada et al. 2020). STM-E exploits all cyclone data drawn from a specific geographical region of interest, provided that certain modelling conditions are not violated by the data. This means in principle that STM-E provides less uncertain estimates of return values than statistical analysis of cyclone data at a 60 single location. To date however, the STM-E methodology has not been directly validated: the objective of the present work is therefore to provide direct validation of return values (in terms of bias and variance characteristics, for return periods T of hundreds of years) from STM-E analysis using sample data for modelling corresponding to a much shorter period T 0 (< T ) of observation, drawn from a full synthetic cyclone database corresponding to a very long period T L (T L > T ) of observation.
In the following sections, we present a motivating application in the region of the Caribbean archipelago of Guadeloupe, for 65 which synthetic cyclone data are available for a period T L corresponding to more that 3,000 years. We use the STM-E method to estimate the T = 500-year return value for SWH, and its uncertainty, based on random samples of tropical cyclones corresponding to T 0 = 200 years of observation. We compare estimates with empirical maxima from random samples corresponding to T years of observation from the full synthetic cyclone data (covering T L years), and from standard extreme value estimates obtained using data (corresponding to T 0 years) from the specific location of interest only. Section 2 provides an outline of the 70 motivating application. Section 3 describes the STM-E methodology. Section 4 presents the results of the application of STM-E to the region of main island pair (Basse-Terre and Grande-Terre) of Guadeloupe. Discussion and conclusions are provided in Section 5. The orange and green tracks respectively represent those of Maria (2017) and Hugo (1989) cyclones (data extracted from Landsea and Franklin (2013) with cyclone status "Hurricane"). Right: Enlarged view of the diagnostic region. The red rectangle indicates the region in the vicinity of Guadeloupe archipelago where the return levels are estimated. within 400km of the study area on average for the period 1970-2019. Almost all events emanated from south east. More than 80% of the events passed close to the northern and eastern coasts of Guadeloupe's main island pair. 85 To assess the cyclone-induced storm surge hazard, Krien et al. (2015) set up a modelling chain similar to that described in the introduction: they randomly generate cyclonic events using the approach of Emanuel et al. (2006), and compute SWH and total water levels for each event over a wide computational domain (45-65W, 9.5-18.3N) using the ADCIRC-SWAN wave-current coupled numerical model. The interested reader can refer to Krien et al. (2015) for more implementation and validation detail.
In the present work, we use a total of 1971 synthetic cyclones passing nearby Guadeloupe (representative of 3,200 years, i.e. 90 about 0.6 cyclone per annum) and the corresponding numerically calculated SWH. These results are used to derive empirically the 100-year SWH around the coast of Guadeloupe's main island pair for a smaller area of interest (60.8 − 62.0 • W, 15.8 − 16.6 • N; see Figure 2). These results are useful to assess flood risk at local scale, since they provide inputs of high resolution hydrodynamic simulations (see e.g. the use of wave over-topping simulations at La Reunion Island by Lecacheux et al. 2021).
In the following, we analyse extreme SWH at 19 coastal locations around Guadeloupe's main island pair (on the 100m iso-95 depth contour, see blue stars in Fig. 2), and at 12 locations along a line transect emanating to the north east from the island, corresponding to increasing water depth (see red triangles in Figure 2).
To illustrate the SWH data, Figure 3 depicts the spatial distributions of maximum SWH per location for the four cyclones with the largest single values of SWH in the whole synthetic cyclone database. All cyclones propagate from the south-east to the north-west with intense storm severity near the cyclone track, which reduces quickly away from the track.

Methodology
In this section we describe the STM-E methodology used to estimate return values in the current work. Section 3.1 motivates the STM-E approach, and Section 3.2 outlines the modelling procedure. Section 3.3 provides a discussion of some of the diagnostic tests performed to ensure that STM-E modelling assumptions are satisfied.

Motivating the STM-E model 105
The STM-E procedure has been described in Wada et al. (2018) and Wada et al. (2020). The approach is intended to provide straightforward estimation of extreme environments over a spatial region, from a relatively small sample of rare events such as cyclones, the effects of any one of which do not typically influence the whole region. For each cyclone event, the spacetime characteristics of the event are summarised using two quantities, the space-time maximum (STM) of the cyclone and the spatial exposure (E) of each location in the region to the event. For any cyclone, the STM is defined as the largest value of SWH  The key modelling assumptions are then that (i) the future characteristics of STM and E over the region will be the same as those of STM and E during the period of observation, (ii) in future, at any location, it is valid to associate any simulated

STM-E procedure
The steps of the modelling procedure are now described. The first three steps of the procedure involve isolation of data for the analysis. (a) An appropriate region of ocean is selected. The characteristics of this region need to be such that the underpinning 120 conditions of the STM-E approach are satisfied (as discussed further in Section 3.3). (b) Then for each tropical cyclone event occurring in the region, the largest value of SWH observed anywhere in the region for the period of the cyclone (STM) is retained. (c) Next, per location in the region, the largest value of SWH observed during the period of the cyclone, expressed as a fraction of STM, is retained as the location exposure E to the cyclone.
The next three steps of the analysis involve statistical modelling and simulation. (d) First, an extreme value model is esti-125 mated using the largest values from the sample of STM; typically, a generalised Pareto distribution (see e.g. Coles et al. 2001) is assumed. Then a model for the distribution of location exposure E is sought; typically we simply re-sample at random with replacement from the values of historical E for the location, although model-fitting is also possible. (e) Next, realisations of random occurrences of STM from (d), each combined with a randomly-sampled exposure E per location, permits estimation of the spatial distribution of SWH corresponding to return periods of arbitrary length. (f) Finally, diagnostic tools are used to 130 confirm the consistency of simulations (e) under the model with historical cyclone characteristics.

Diagnostics for STM-E modelling assumptions
The success of the current approach relies critically on our ability to show that simplifying assumptions regarding the characteristics of STM and exposure are justified for the data to hand. In particular, the approach assumes that (i) the distribution STM does not depend on cyclone track, environmental covariates, space and time, and (ii) the distribution of exposure per location 135 does not depend on STM, cyclone track, environmental covariates and time. Diagnostic tests are undertaken to examine the plausibility of these conditions for region of ocean of interest for each application undertaken. Establishing the validity of the STM-E conditions is vital for credible estimation of return value estimates. Section 5 of Wada et al. (2018) provides a detailed discussion of some of the diagnostic tests that should be considered to judge that the STM-E conditions are non violated in any particular application.

140
The absence of spatial trend in STM over the region can be assessed by estimating the linear or rank correlation between STM values along longitude-latitude transects with arbitrary orientation in the region. For comparison, the "null" distribution of correlation can be estimated using random permutations of the STM values.
Return value estimates from STM-E are also potentially sensitive to the choice of region for analysis. We assume that the extremal behaviour of STM can be considered homogeneous in the region, suggesting that the region should sufficiently 145 small that the same physics is active throughout it. However, the region also needs to contain sufficient evidence for cyclone events and their characteristics to allow reasonable estimation of tails of distributions for SWH per location. The absence of dependence between STM and E per location can be assessed by calculating the rank correlation between STM S (a space-time maximum) and exposure E j (at location j) using Kendall's tau statistic. If the values of S and E j increase together, the value of Kendall's tau statistic will be near to unity. If there is no particular relationship between S and E j , the value of Kendall's tau

Modelling STM and estimating return values
Suppose we have isolated a set of n 0 values of STM using the procedure above. We use the largest n ≤ n 0 values {s i } n i=1 , corresponding to exceedances of threshold ψ n , to estimate a generalised Pareto model for STM, with probability density with shape parameter ξ ∈ R and scale parameter σ n > 0. Choice of n is important, to ensure reasonable model fit and biasvariance trade-off. The estimated value of ξ should be approximately constant as a function of n for sufficiently small n. The 165 full distribution F S (s) of STM can then be estimated using where F * n (s) is an empirical "counting" estimate below threshold ψ n , and τ n is the non-exceedance probability corresponding the ψ n , again estimated empirically.
Using this model, we can simulate future values H j of SWH at any location j, (j = 1, 2, ..., p) in the region, relatively 170 straightforwardly. Suppose that E j is the location exposure at location j, and F Ej its cumulative distribution function, estimated empirically. Since then H j = E j × S, the cumulative distribution function of H j can be estimated using where f S (s) is the probability density function of STM, corresponding to cumulative distribution function F S (s) estimated in Equation 2.

Application of STM-E to cyclones SWH near Guadeloupe
The STM-E methodology outlined in Section 3 is applied to data for the neighbourhood of Guadeloupe's main island pair

Details of STM-E application
The spatial region of interest is the neighbourhood of Guadeloupe's main island pair in the Caribbean Sea, corresponding to 190 approximately longitudes 12 • − 18 • N and latitudes 58 • − 65 • W (see Figure 1). An initial analysis using Kendall's tau suggests shows dependency of STM and exposure, with stronger cyclones tending to pass through the western part of the region. However, if a very high threshold ψ ≈ 20m was selected for analysis, reasonable decoupling of STM and E could be achieved, with relatively less intense tropical cyclones neglected. Since the focus of the current work is the ocean environment of Guadeloupe archipelago, a smaller region (see Figure 1, right panel) was defined.

195
For this region, Kendall's tau indicated low dependence between STM and E for thresholds ψ of 10m and above, as illustrated in the left panel of Figure 4.
The right panel of Figure 4 shows the location and magnitude of STM for each of the n = 60 largest cyclones observed in the region. There is no obviously no spatial dependence between the size of STM and its location. In Section 3.3, we discuss the use of rank correlation of STM along latitude-longitude transects as a means to quantify dependence in general. In fact, the  east of the main island pair of Guadeloupe. We focus on the north-east part, because it is where the exposure to cyclones is the highest. The contour and transect are illustrated in Figure 2, and location numbers are listed.
Focus of the analysis is estimation of the T = 500-year return value for SWH on the iso-depth contour and line transect, based on T 0 = 200 years of data, to quantify the uncertainty in the 500-year return value using STM-E analysis, the follow-215 ing procedure based on random sub-sampling was adopted. For STM-E analysis, the following procedure was adopted.

Benchmarking against the full cyclone database and single-location analysis
One important advantage of the current synthetic cyclone database is that it corresponds to long time period relative of T F = 3, 200 years, much longer than the return period of T = 500 years being estimated in the current analysis. Thus, we are able to estimate the 500-year return value at any location using the full synthetic cyclone data, by simply sampling 310 cyclones 230 (corresponding to 500 years, with annual cyclone occurrence rate of 0.62) at random. This provides a direct empirical estimate.
From previous work, a key advantage found using the STM-E approach is that it provides less uncertain estimates at a location compared with conventional "single location" analysis. We wish to demonstrate in the current work that this is also the case. For this reason, we also calculate estimates for comparison with those from STM-E, based on independent analysis of cyclone data from each location of interest. The procedure is as follows. single-location (red) and full synthetic cyclone data (black); in general, there is good agreement between estimates per location.
It can be seen however (from the longer red whiskers) that the uncertainty in single-location estimates is greater than estimated from the full data and from STM-E. In general, the uncertainty in return value estimated using STM-E more closely reflects that seen in the full synthetic cyclone data.
Per location, as the number points used for STM-E estimation increases, there is evidence for reduction in the uncertainty 250 with which the return value is estimated, as might be expected. However, there is also some evidence for a small increase in the mean estimated return value. This is explored further in Section 4.5. There is no corresponding evidence for reduced uncertainty in the single-location analysis. There are more outlying estimates of return value for single-location analysis than for STM-E. Whiskers represent the 2.5% to 97.5% interval. Exceedances of this 95% interval are shown as dots. The green dashed line and the right-hand y-axis give the water depth at each location.
The corresponding results for the line transect analysis using maximum likelihood estimation is given in Figure 7. The 255 general characteristics of this figure are similar to those of Figure 6. The return value increases as would be expected with increasing water depth. Single-point estimates are more variable that those from STM-E. There is again some evidence that the STM-E median estimate increases with increasing sample size, but this trend is small compared to the uncertainty in the return value estimated empirically from the full synthetic cyclone data set. We infer from the analysis that water depth has little effect on the performance of the STM-E approach.

Results estimated using probability weighted moments
Estimates for the 500-year return value on the iso-depth contour, obtained using the method of probability weighted moments, are shown in Figure 8. The behaviour of STM-E and single-location estimates shown is very similar to that illustrate for maximum likelihood estimation in Figure 6. There is some evidence that both STM-E and single-location estimates obtained using probability weighted moments show somewhat lower bias. Again, this is perhaps not surprising given the relative characteris-265 tics of maximum likelihood and probability weighted moment estimators for small samples. This is further discussed in Section 4.5.
Results for the line transect using probability weighted moments are given in Figure 9. Again, the figure shows similar trends to Figure 7, except perhaps for a slightly lower bias in estimates using probability weighted moments.

Assessment of model performance 270
Results from Sections 4.3 and 4.4 suggest that there is a small increasing trend in return value estimates from STM-E as a function of increasing sample size for inference. Results also suggest that the method of probability weighted moments gives less bias in estimates for return values. We investigate these effects further here. Figure 10 gives estimates for the 500-year return value of space-time maximum STM (as opposed to the full STM-E estimate) as a function of sample size used for estimation, using maximum likelihood estimation (upper panel) and probability weighted moments. Also shown (in red, at 275 zero on the x-axis) is the empirical estimate of the 500-year STM return value obtained directly from the synthetic cyclone data. The figure shows a number of interesting effects. Firstly, maximum likelihood estimates are in this case biased low, and the extent of the bias reduces with increasing sample size n. It should be noted, however, that the extent of the bias is small relative to the uncertainty in the the return value. Secondly, the method of probability weighted moments is less biased for small sample sizes. But interestingly, the mean (or median) return value estimate also appears to increase with sample size, so that at 280 a sample size of 60, probability weighted moment estimates are clearly biased high in this example. Further, uncertainties on return value estimates from probability weighted moments are larger than those from maximum likelihood estimation.
There are many studies in the literature comparing the performance of different methods of estimation of extreme value models. The method of probability weighted moments is known to perform relatively well relative to maximum likelihood estimation for small samples (see, e.g., Jonathan et al. (2021), Section 7 for a discussion). For small samples, for example, Figure 9. 500-year return value for SWH estimated using probability weighted moments for the line transect. The x-axis gives the reference numbers of the 12 locations on the transect. Briefly, for each location, blue and red box-whiskers summarise the estimated return value from STM-E and single-location analysis respectively; see caption of Figure 6 for other details. The black box-whisker plot per location corresponds to the empirical estimate of the return value obtained directly from the synthetic cyclone data. The green dashed line and the right-hand y-axis give the water depth at each location. maximum likelihood estimation is known to underestimate the generalised Pareto shape parameter, and over-estimate the corresponding scale parameter, leading to bias in return value estimates. Given these facts, we do not consider the small trends shown in the current analysis to be particularly interesting, or of particular concern. However, the trends do illustrate the importance of performing an extreme value analysis with great care, particularly for small sample sizes.
One of the assumptions underpinning the STM-E approach is that the exposure distribution at a location is not dependent on 290 the magnitude of STM. We investigate this effect in Figure 11. The figure gives exposure curves for each of the 100 random data samples of 200 years used to estimate the STM-E model and return value characteristics in this section, for location 10 (see Figure 2). Panels correspond to different sample sizes. We are concerned in particular that there are no patterns in exposure distribution related to the size of STM for the corresponding data sample. For this reason, each panel also shows the exposure distribution corresponding to the data sample with largest and smallest STM in red and blue. It can be seen that the coloured 295 curves are not obviously unusual with respect to the other curves. Similar plots for other locations also do not give cause for concern.

Discussion
This work considers the estimation of P -year return values for SWH over a geographic region, from small sets of T 0 years of synthetic tropical cyclone data, using the STM-E (space-time maxima and exposure) methodology. We assess the methodology   For reasonable application of the STM-E approach, it is important that characteristics of tropical cyclones over the region under consideration satisfy a number of simplifying conditions. These conditions are shown not to be violated for a region around Guadeloupe archipelago, but that use of the STM-E method over a larger spatial domain would not be valid (see e.g. Wada et al. 2019). This demonstrates that selection of an appropriate geographical region for STM-E analysis is critical to its 310 success. Once such a region is specified, we find that the space-time maxima and exposure (STM-E) method provides a simple but principled approach to return value estimation within the region from small samples of tropical cyclone data. Return value estimates from STM (see e.g. Figure 10) show a small increasing bias with increasing sample size n for extreme value estimation. However, the resulting bias in full STM-E return values is small, especially compared with the inherent variability of the return value estimate. Corresponding estimates based on single-location analysis also show relatively small 315 but increasing negative bias with increasing sample size. In the present work, the tail of the distribution of STM was estimated by fitting a generalised Pareto model, using either maximum likelihood estimation or the method of probability weighted moments. Estimates for extreme quantiles of STM using either approach are in good agreement. Table 1 summarises the performance of STM-E and single-location analysis in estimation of the bias and uncertainty of the 500-year return value, relative to empirical estimates from the full synthetic cyclone data, for analysis sample sizes of n=20, 30, Here,h( ; n, Mth) andh 0 (n) correspond to the median 500-year return value estimated using sample size n from inference method Mth (either maximum likelihood or probability-weighted moments) and directly from the full synthetic cyclone STM-E is less biased than single location estimates except for sample sizes 20 and 30 using probability weighted moments.
Bias from STM-E is smaller in magnitude by more than 1.0m for sample size n ≥ 50 using probability weighted moments. In 330 terms of uncertainty, STM-E always produces more realistic estimates of return value uncertainty that single location analysis.
For sample sizes n = 50 and 60, STM-E only over-predicts the inherent uncertainty of the return value by some 10-30%. In conclusion, if a single location approach had been adopted for the current application, the resulting uncertainty over-estimation would have been more than 100%, and the bias could have reach values up to 1.50m. We hope that the table provides some general guidance about the expected future performance of STM-E for locations similar to Guadeloupe archipelago, compared 335 with single location analysis.
It is important to keep in mind that an appropriate choice of sample size n for STM-E analysis is likely to be related to the size n 0 of the full sample available, and the period T 0 to which the sample corresponds. For example, in the current work, n = 20 (and n 0 = 124) is approximately equivalent to the largest 15% of cyclones for the sample period T 0 = 200 year. That is, the smallest cyclone considered in the n = 20 STM-E model has a return period of the order of 10 years. With n = 60, we 340 use approximately half the sample for STM-E analysis, and the smallest cyclone in the STM-E analysis has a return period of the order of 3 years. In another application where T 0 = 50 (and n 0 = 32) say, it is likely that a greater proportion of the sample would need to be used to estimate the STM-E model reasonably.
The table, combined with the validation experiments reported above, confirm the principle intuited by previous works' results (Wada et al. 2018, Wada et al. 2020 Table 1. Performance of STM-E and single-location analysis in estimation of the bias and uncertainty of the 500-year return value, relative to empirical estimates from the full synthetic cyclone data, for analysis sample sizes of 20, 30, 40, 50 and 60. Bias B is assessed as the average difference (over the iso-depth contour and line transect analyses) between the median STM-E (or single-location) estimate and the median estimate from the full cyclone data, for that sample size. Similarly, uncertainty U is assessed in terms of the average ratio of the width of the 95% uncertainty band of the STM-E (or single-location) estimate to the width of corresponding uncertainty band from the full cyclone data, minus one. Thus, zero values of B and U indicate agreement between the model-based inference and that from resampling the full synthetic cyclone database.
ocean (McInnes et al. 2014) or Indian Ocean basin (Lecacheux et al. 2012) where cyclone-induced storm wave data is limited, STM-E achieves low bias and realistic levels of uncertainty, and should be preferred to the single location approach. It confirms that STM-E can be considered a valuable tool for any coastal risk practitioner.