## Abstract

Long-term estimation of extreme wave height remains a key challenge because of the short duration of available wave data, and also because of the possible impact of climate variability on ocean waves. Here, we analyse storm-based statistics to obtain estimates of extreme wave height at locations in the northeast Atlantic and North Sea using the NORA10 wave hindcast (1958–2011), and use a 5 year sliding window to examine temporal variability. The decadal variability is correlated to the North Atlantic oscillation and other atmospheric modes, using a six-term predictor model incorporating the climate indices and their Hilbert transforms. This allows reconstruction of the historic extreme climate back to 1661, using a combination of known and proxy climate indices. Significant decadal variability primarily driven by the North Atlantic oscillation is observed, and this should be considered for the long-term survivability of offshore structures and marine renewable energy devices. The analysis on wave climate reconstruction reveals that the variation of the mean, 99th percentile and extreme wave climates over decadal time scales for locations close to the dominant storm tracks in the open North Atlantic are comparable, whereas the wave climates for the rest of the locations including the North Sea are rather different.

## 1. Introduction

Extreme wave climate is important for the design of offshore structures, marine operations and coastal defences. Accurate estimation of such extremes requires long-term, good-quality measurements, which are rarely available, and hence this remains a key issue. Studies of the extreme wave climate have increased in number in recent years. Both global and regional estimates of extreme significant wave height (*H*_{s}) values, as well as long-term trend analysis, have been produced [1–5]. Seasonal and interannual variability of extreme *H*_{s} values has been assessed, and associated with large-scale pressure anomalies or teleconnections, such as the North Atlantic oscillation (NAO) in the northeast Atlantic [6–9].

Are the properties of waves in the Atlantic over the last 50 years a reliable guide to what may happen in the next 50 or 100 years? We aim to tackle this question by seeking correlations between the variability of the extreme wave heights in the northeast Atlantic and the North Sea and large-scale atmospheric teleconnections. For this region, we include the NAO, and the next two modes, which are the East Atlantic (EA) pattern, and the Scandinavian pattern (SCA). The NAO measures variation in the Atlantic eddy-driven jet that controls the near-surface westerly winds [10], the EA pattern affects the position of the dominant North Atlantic storm track and jet streams [10] and the SCA pattern influences cyclone activity across Northern Europe and Eurasia [11]. Hence, a link between pressure anomalies and extreme wave climate is likely. The NAO and other modes are found to correlate well with mean wave climate [12], so we extend this approach to extreme waves. The temporal pattern of each mode is characterized by a climate index, values of which are obtainable from National Oceanic and Atmospheric Administration (NOAA) Climate Prediction Center (CPC; www.cpc.ncep.noaa.gov), and are based on the rotated empirical orthogonal function analysis by Barnston & Livezey [13]. The available indices cover January 1950 to the present.

Here we present extreme value analysis on storm severity by using a peaks-over-threshold (POT) technique to fit the upper tail of the data with an exponential distribution. The analysis is performed on a sliding window over 5 year blocks to produce temporal variability for correlation with the NAO and other modes. The results are confirmed by comparison with the decadal variation of the 99th percentile of *H*_{s} which requires no extrapolation. This study is a continuation of previous work by Taylor and co-workers [14,15] using earlier hindcast wave models, and Santo *et al*. [12], who looked at annual mean wave power variability using the same hindcast wave model.

## 2. Data and methods

We examine hindcast wave data at 12 locations in the northeast Atlantic and the North Sea from 1958 to 2011, as shown in figure 1*a*. The hindcast data are taken from the Norwegian 10 km Reanalysis Archive (NORA10), for more information, see [16]. We also use measured buoy data for Haltenbanken from 1980–2000 (and Forties from 1974 to 1996 but this contains large gaps, so analysis of decadal variation of extreme values is not attempted). Both hindcast and measured wave data were sampled at 3 h intervals, and contain information such as date, time, significant wave height (*H*_{s}), peak spectral wave period (*T*_{p}), wind speed, wind and wave directions. For the extreme wave height analysis, *H*_{s} and *T*_{p} are used. Previous analysis of mean wave power variability using the same data had investigated the comparisons of hindcast/buoys in terms of the annual mean values at Haltenbanken and Forties, and in general, the agreement is reasonable [12].

Figure 1*b* shows directional wave rose distributions for the largest 1000 storms at three locations: St Kilda, Bruce and Valhall (points 4, 8 and 12, respectively), shown in terms of the most probable maximum individual wave height (*H*_{mp}), which is a measure of storm severity (discussed in §2a). The wave direction for each *H*_{mp} is set to the wave direction of the peak *H*_{s} within each storm. The range circles show the frequency (in %) of *H*_{mp} occurrence normalized by the total occurrences for all directional bins of 10°. The wave rose distributions for the rest of the locations are provided in electronic supplementary material, figure S1.

The dominant large waves in the northeast Atlantic locations (points 1–7, including St Kilda) come from the west generated by storms in the open North Atlantic, and consistent with the dominant wave direction contributing to the mean wave power variability, see [12]. For the North Sea locations (points 8–11) except Valhall (point 12), the dominant large waves come from effectively two directional sectors (northwest and southeast), similar to the two dominant wave directions contributing to the mean wave power variability. It is also interesting to see that, at Bruce, waves propagate from the west through the Orkney–Shetland gap and from the north coming down the east coast of the Shetlands.

Northwesterly waves get weaker further south. Thus, one would expect a systematic decrease in the strength of the extreme waves moving south down the North Sea. There is a second contribution of waves from effectively the opposite direction (from the southeast), justifying partitioning the data on incoming wave direction. Curiously, Valhall is the only location in the North Sea that we have examined where the dominant large waves are from one directional sector (from north to southwest); there are smaller waves from the southeast sector but these are not shown as they are all below the POT-threshold.

### (a) Storm-based identification

Simple extreme value analysis methods require observations to be independent and identically distributed. Successive *H*_{s} values are strongly correlated as they are likely to be part of a single longer storm. Hence, storms are considered for long-term extrapolation instead. The entire (*H*_{s},*T*_{p}) record is grouped into storms, which are effectively uncorrelated. Only winter storms are considered (from October to March every year). A single storm at a fixed location is assumed to last not more than 24 h, a reflection of the time scale for the motion of winter depressions in the northeast Atlantic and the North Sea [17]. Because the probability distribution of the upper tail of the individual wave height distribution for a storm is governed solely by the portion which has *H*_{s}>0.8*H*_{sMax} [18], we use this threshold to define a single storm. We investigate the sensitivity of this threshold and the assumed upper limit for storm duration of 24 h on the 100 year return period estimates of *H*_{mp}, and confirm that the results are robust, see electronic supplementary material, table S1 and figure S2.

Once all storms have been identified, the most probable largest individual wave in each storm (*H*_{mp}) could be obtained following the method by Tromans & Vanderschuren [17], using convolution integrals. A simpler process to obtain *H*_{mp} by random sampling is presented as follows. For each (*H*_{s},*T*_{p}) record in a storm, the corresponding zero crossing period of individual waves, *T*_{z}, is estimated using *T*_{z}=0.779*T*_{p} assuming the JONSWAP spectral shape, which gives an estimate of the number of waves (*N*_{w}) in each 3 h sea-state [19,18]. Assuming a Rayleigh distribution to account for the intrasea-state and intrastorm variability of individual wave heights within a sea-state, the probability (*P*) that all *N*_{w} peaks are simultaneously less than a wave amplitude, *A*, can be expressed as
*σ*=*H*_{s}/4 and *A* is obtained by randomly sampling the probability *P* between 0 and 1. This is performed for each *H*_{s} record (each 3 h sea-state) in a storm and *σ*_{s} and *N*_{s}) performed by least-squares minimization. These parameters are those of an equivalent rectangular area storm, with *σ*_{s} being the intensity and *N*_{s} being the measure of storm duration. Other than using the fit to identify *H*_{mp}, we make no other use of them here.

The peak value of the empirical PDF (or the mode) is the most probable maximum amplitude of the wave (*A*_{mp}). For linear wave theory, the most probable maximum wave height within a storm is *H*_{mp}=2×*A*_{mp}. The use of *H*_{mp} accounts for the possibility that the largest wave may not always occur in the most severe 3 h sea-state and storms often last for more than 3 h. Hence, it is a more robust measure of storm severity than *H*_{sMax} [14,15,17]. While *H*_{mp} is a convenient parameter for estimating storm severity, in defining it, we are making a narrow-banded assumption about the individual waves in each sea-state. In addition, the 100 year return period estimate of individual wave height (*H*_{mp} should not be used as a design wave height. Nevertheless, because *H*_{mp} is less noisy than

As an alternative to the Rayleigh distribution, we could have used the distribution for wave height proposed by Forristall [20]. This has the drawback that wave height as conventionally defined is not a point-wise property of the record (being the sum of a crest height and following trough in time, with these separated by approximately half the wave period). Forristall's proposed distribution for the short-term statistics of wave crests [21] avoids this problem, but does not include the second-order difference term. In addition, crest elevation is a combination of a linear wave amplitude and a mostly double frequency sum contribution that makes crests and troughs different shapes. For structural reliability calculations, engineers usually start from the linear wave amplitude, or twice this being the local height of the wave envelope as used here; these being satisfactorily (and slightly conservatively) approximated using the simple Rayleigh distribution.

While we regard *H*_{mp} as an important storm-based parameter for engineering applications, we also present the results of extreme value analysis on the more commonly used *H*_{sMax} by using the same storm-based identification and taking the largest *H*_{s} record within each storm as *H*_{sMax} (one entry per storm). It is observed that *H*_{mp}∼2×*H*_{sMax} for 100 and 1000 year return period estimates. We also present the results on the 99th percentile of the complete *H*_{s} record over all winters (from October to March every year) to produce an in-sample estimate for comparison with *H*_{mp} and *H*_{sMax}. Unless otherwise specified, most of the results are presented in terms of 100 year return period of *H*_{mp}, but we will show that the other parameters have very similar behaviour over decadal time scales.

### (b) Peaks-over-threshold technique

The identified storms are ordered in terms of *H*_{mp} (and *H*_{sMax}) magnitude, and the POT technique is used for the long-term extrapolation. We seek a tail form for the cumulative density function (CDF) that fits the observations above a threshold (or exceedances). Formally, this is described by the generalized Pareto distribution (GPD) [22]. In general, there are three simple characteristic forms of the GPD described in terms of the thickness of the high tail: ‘fat’ tail, ‘thin’ tail and tail with an upper limit.

Here we assume that the exceedances fall into the class of ‘thin’ asymptotic tails with an exponential form (see eqn. 4.4 in [22]), with the CDF written as
*H*_{c} is the cut-off (threshold) value of *H*_{mp}, and *α* is the scale parameter, to be obtained by the maximum-likelihood method.

We considered the possibility of a data transformation prior to estimating the exceedances (as discussed by [23]). This would be equivalent to replacing the exponential form with a Weibull fit commonly used by oceanographers (see eqn 6.8-4 in [18]), but for this particular case, the simple exponential fits well. Arguably, the exponential fit could be assumed across all the locations, because most of the storm-driven large waves are coming across the North Atlantic (from the wave rose distributions). It should be noted that most of the analysis is aimed at the 100 year (10^{−2} yr^{−1}) wave. With 54 years of hindcast data from NORA10, little extrapolation is required. Therefore, we expect that almost any type of reasonable fit would yield comparable results.

## 3. Extreme wave height

First, we perform extreme value analysis using the entire 54 year hindcast storm record from 1958 to 2011, which assumes a stationary wave climate over 54 years, ignoring any long-term trends. Subsequently, to account for non-stationarity, we apply a sliding window over 5 year blocks assuming the data are stationary and identically distributed only within each 5 year period, and produce estimates of the temporal variability of the 100 year *H*_{mp} (10^{−2} risk per year) for correlation with the teleconnections. The stationarity assumption for the 54 year record is perhaps reasonable given that the long-term trend from the sliding window analysis is small, despite the large decadal variability that will be apparent later. It is worth stressing that stationarity is commonly assumed in standard offshore industry practice, unless intrarecord variability is specifically looked for. Using a POT method, the twin requirements of enough large values within the data and the sensitivity of the extrapolated values to threshold, so extreme value theory can be assumed, both drive the inclusion of as much data into a single analysis as possible.

Because robust threshold behaviour is important for POT analysis, we investigate the sensitivity to threshold level using the entire 54 years of record. We use the number of storms (*N*) as a threshold for every location, hence the value of *H*_{c} will vary for each location depending on the storm severity over time. Figure 2*a* shows the variation in the predicted 100 year *H*_{mp} value at Schiehallion as the threshold *N* (and *H*_{c}) for the fit is changed. The value of the 100 year *H*_{mp} is reasonably stable around *N*∼1000 largest storms, and this generally holds elsewhere, suggesting the choice of threshold level is robust. A threshold level of *N*∼1000 storms corresponds to approximately 20 storms per winter, roughly one to two per week.

Figure 2*b* shows the long-term fit of the *H*_{mp} for Schiehallion. Every data point represents an individual storm severity (*H*_{mp}) from the hindcast, and all these lie within or close to the envelope, or trumpet (estimated 5–95% confidence interval), which is obtained by a parametric bootstrap. A very similar trumpet is also obtainable by non-parametric bootstrap (by sampling from the data itself). The mean of the bootstrap fits very closely to the original data fit, and the resultant trumpet is close to symmetric around the mean value, suggesting that the analysis is robust. Generally, the variation of the width of the trumpet at the 5–95% level for the 100 year *H*_{mp} is ∼7% of the estimated mean value.

The *H*_{mp} fits for the rest of the locations are summarized in table 1 and plotted in electronic supplementary material, figures S3 and S4. The ratio of the 1000 to 100 year *H*_{mp} (a measure of the slope of the long-term extrapolation) is about 1.17–1.18, consistent for all locations. In general, for the open North Atlantic locations, St Kilda has the most severe extreme wave height values, followed by Corrib, Schiehallion and Haltenbanken. For the North Sea locations, the extreme wave heights decrease slightly from north to south, consistent with the largest *H*_{sMax} values in the hindcast and with the wave height rose distributions in §2.

Figure 3*a* shows the 5 year sliding window of the 100 year *H*_{mp} for Schiehallion. The definition of the storm year is based on the year of December of the third winter, and two winters either side. With *N*∼1000 storms over 54 years, a 5 year sliding block contains ∼100 storms (on average), which we judge sufficient for extreme value analysis. The signal is dominated by significant decadal variability with no obvious long-term trend, and it is this decadal variability that will be correlated with the teleconnections. Also shown is the long-term mean value (dashed line) of the 5 year based 100 year value, which here for Schiehallion is 34.08 m, only 2% different from the value for the entire 54 years storm record (34.73 m, table 1). Thus, the long-term trend of the data is small, and this also holds for all locations. The degree of the temporal variability for all location ranges from 5% to 8% in terms of coefficient of variation.

On the same figure, the estimated 5–95% confidence interval represents the variability by random sampling within each 5 year block. Given that the long-term mean value is within the broad 5–95% confidence interval, is the time-varying signal driven by real large-scale atmospheric variability, or is it just statistical noise from a long-term constant distribution? In this paper, we present evidence that this variation is real because of the strong correlation with the NAO and other modes.

Figure 3*b* shows the 5 year sliding window of the 100 year *H*_{sMax} (i) and the 99th percentile of *H*_{s} (ii) for Schiehallion. The sliding window results for the 100 year *H*_{sMax} is obtained using the same procedure as for the 100 year *H*_{mp}, whereas for the 99th percentile of *H*_{s}, the value of the 99th percentile is obtained on an annual basis over 54 years, and a 5 year moving average is run across the signal to produce comparable 5 year sliding window results. In general, the decadal variability in both the 100 year *H*_{sMax} and the 99th percentile of *H*_{s} is similar to that of the 100 year *H*_{mp}. However, there is some slight difference in the year-to-year variation of the 99th percentile of *H*_{s} when compared with the other two results, which is likely owing to finite size sample variability in the POT analysis. This is also evident from the width of the estimated 5–95% confidence interval which is narrower for the in-sample estimate of the 99th percentile of *H*_{s}. In addition, it is worth noting that *H*_{sMax}∼2× 99th percentile of *H*_{s}, so we are looking at a different part of the distributions which are treated in a rather different manner, hence the level of agreement in the ‘wiggles’ is encouraging. We consider that the similarity in the decadal variability of the POT results with that of the in-sample estimate demonstrates the robustness of our extreme value analysis.

Hindcast/buoy comparison in terms of individual *H*_{s} above 5 m is performed for Haltenbanken and Forties. Normalized root-mean-squared error is used to quantify the comparison: 15% for Haltenbanken and 11% for Forties with no bias. For storm-by-storm analysis, the comparison is performed for Haltenbanken for 21 years, from 1980 to 2000, with *N*∼450 storms to yield a consistent ∼100 storms per 5 year block on average. We compare the largest 450 *H*_{sMax} and *H*_{mp} values (ranked order storm), and the agreement is satisfactory (figure 4). The comparison of the 100 year and 1000 year extreme values are close (table 1). For a 5 year sliding window on the same period of 21 years, the temporal variation of the 100 year values are similarly close, with *R*^{2}=0.75 (figure 5, buoy-based results in red), supporting the use of the hindcast wave data.

## 4. Correlation with the teleconnections

Following the method of correlating the annual mean wave power with the teleconnections described in [12], we correlate the 5 year sliding window estimates of the 100 year *H*_{mp} with the NAO, the EA and the SCA using the same winter average of the climate indices (obtained by averaging six month values of monthly teleconnection indices centred around the middle of January) and the same form of the predictor model, which is based on linear regression. We then compare the performance of this three-term predictor model with an extended six-term predictor model, which incorporates both the three climate indices and their Hilbert transforms of the indices (which represent 90° phase shifted versions and are related to the time derivative of the indices without amplification of high-frequency noise) [24,25].

The form of the six-term predictor model is expressed as
*H*_{mp} over the period of available hindcast data, *NAO*(*t*) is the temporal winter average of the NAO index, *NAO*(*t*), likewise for *EA*(*t*) and *SCA*(*t*). The subscript *H* represents the Hilbert transform of the indices. *b*, *c* and *d* are non-dimensionalized constants resulting from the variance minimization; these reflect the relative importance of the NAO, the EA and the SCA signals (subscripts 1), and their Hilbert transform counterparts (subscripts 2) in describing variability of extreme wave height. For *b*_{2}=*c*_{2}=*d*_{2}=0, the above form reduces to the three-term predictor model. The predictor model is trained over the period of available hindcast data.

The Hilbert transform is performed on each index running from 1950 to present, which has been patched at each end of the signal with a scaled autocorrelation of the index to minimize end effects. Arguably, we are not over-fitting the model as the additional three index signals are derived analytically from the original three indices (a Hilbert transform of a signal is orthogonal to the original signal). For both predictor models, we impose a low-pass filter of 5 years to all climate indices and Hilbert-transformed indices (moving average), to be consistent with our 5 year sliding window analysis.

The correlation results of the 5 year sliding window for all the locations are summarized in table 1. In general, the correlations at all locations are improved with the six-term predictor model compared with the three-term model, demonstrating the ability of the six-term model in capturing more of the short-term variability within the training period (50 years). More importantly, the improved results are more consistent, i.e. we obtain reasonable correlations with the NAO and other modes uniformly across the locations in the open North Atlantic with the six-term model. Additional information regarding the optimized non-dimensionalized constants is provided in electronic supplementary material, tables S2 and S3. In general, the NAO has the largest contribution when the *R*^{2} is significant (locations 1–4), with contributions being measured in terms of the relative strength of the non-dimensionalized constants in each predictor model. However, some of the correlation results look peculiar, such as from locations 5–7 (going southwards) where the SCA is seen to have the most dominant contribution over the NAO and the EA.

Figure 5*a* illustrates, using the six-term predictor model, the correlation result of the 5 year sliding window for St Kilda, Haltenbanken and Valhall in terms of the one in 100 year *H*_{mp} (*R*^{2}=0.54−0.67). We note the similarity of the extreme *H*_{mp} and the predictor model over time scales longer than 5 years. This becomes evident when a 5 year moving average is run across each signal without doing any further computation, and the correlation is improved (*R*^{2}=0.69−0.87), see figure 5*b* and the second column under each predictor model in table 1. It is likely that the long-term variability beyond 5 years can be represented well by the NAO and other modes, but for intradecadal variability over less than 5 years, finite size sample variability in the POT analysis introduces significant numerical noise. The correlation results for the rest of the locations are provided in electronic supplementary material, figures S5 and S6 for the six-term model. Analysis of the residual term for each comparison reveals no significant autocorrelation structure beyond a 5 year lag. The correlation result of the one in 100 year *H*_{sMax} for all three locations is similar to that of the *H*_{mp} for the same return period, as evident from the similarity in the temporal variation of the 5 year sliding window result of both as shown previously in figure 3. The correlation result of the 99th percentile of *H*_{s} is improved for Haltenbanken and St Kilda with *R*^{2}=0.83 and 0.71, respectively, whereas the correlation for Valhall is similar to that of *H*_{mp} with *R*^{2}=0.65.

Figure 5*c* shows the variation of the 5 year moving average of the three atmospheric indices used in the predictor model (the associated Hilbert transforms are not shown). One can observe close resemblance of the variation between the NAO and the one in 100 year *H*_{mp} at the locations considered in figure 5*a*, likewise for the variation in the EA, whereas it is less obvious for the SCA. When the NAO is in positive phase (or less negative), the extreme wave activity is observed to increase accordingly, such as for the periods between 1970–1980 and in the 1990s, and vice versa when the NAO is in negative phase (or less positive). The strong relationship between the phase of the NAO and the wave variability, consistent with [26], is supported by the general observation that the positive phase of the NAO (and in combination with other atmospheric modes) is associated with warm and wet winters (plenty of large waves) in Northern Europe, while the negative phase is associated with cold and dry winters (fewer large waves).

As suggested by one of the reviewers, the question of whether the predominantly North Sea predictors EA and the SCA would benefit reconstruction of waves in the North Atlantic, we investigate the relative importance by dropping out the least significant atmospheric mode which is either the EA or the SCA. The effect on the *R*^{2} is assessed, which is shown in table 1 in the third column for each predictor model. The smallest contribution is from the EA for almost all locations, except locations 2–3 where the SCA has the least contribution. In general, the reduction in the *R*^{2} (or Δ*R*^{2}) is larger for the six-term model compared with the three-term model for all locations. From the six-term model, we observe that all three atmospheric indices are significant (Δ*R*^{2} = 5–67%), so should be included in the predictor model except for locations 2–4 (Δ*R*^{2}<5%). Hence, we retain all three indices (and the associated Hilbert transforms) for the subsequent analysis.

From table 1, reasonable correlation (*R*^{2}>0.6) on the extreme wave height is obtained for the open North Atlantic locations (points 1–6) and Valhall (point 12), but the correlations for the rest of the locations (point 7–11) are rather weak. The correlation for Bruce and other locations in the North Sea is weak owing to the existence of two (or more) directional sectors. We split the omnidirectional waves into two directional sectors based on incoming wave directions, and investigate the temporal behaviour of each sector using the same methodology. The waves at Bruce (point 8), for instance, are split into the northwest (230°–50°) and southeast (50–230°) sectors, where 0° is pointing to the north and 90° to the east. Dramatically improved correlation in terms of 5 year moving average is obtained for the northwest sector (*R*^{2}=0.90) when compared with the omnidirectional signal (*R*^{2}=0.35); however, the southeast sector remains mostly unexplained by the NAO and other modes (*R*^{2}=0.41). The correlation results for locations in the North Sea are summarized in table 2, and plotted in electronic supplementary material, figure S7. The correlation for Cornwall is weak perhaps because of the close proximity to shore and long distance sheltering effects from Ireland.

Overall, with the six-term predictor model, we observe that the westerly extreme wave heights at all locations are significantly correlated to the NAO and other modes (*R*^{2}=0.6–0.9). In particular, points 1–4 and 8 have the strongest correlation with the NAO, apart from point 12. Interestingly, these locations are in close proximity to the dominant storm tracks. Because the NAO is known to be correlated with movements of the storm tracks, and latitude shifts of storm tracks are associated with changes of the largest waves in the North Atlantic [27], a strong relationship between the NAO and the extreme wave climate at these locations is expected.

## 5. Reconstructions of extreme wave climates

Significant correlation between the extreme wave climate and the dominant teleconnections over the period of the available hindcast data allows reconstruction of historic extreme wave climate, under the assumption that such relations hold into the past. We use the same proxy climate indices from 1659 to 1998 as used previously in [12]. The proxy indices were obtained by regressing the known indices from the CPC with the historical reconstructed monthly 500 mbar pressure maps computed by Luterbacher *et al*. [28], over the overlapping period from 1950 to 1998, see [12] for more details. Because our reconstruction analysis relies solely on the pressure data reconstructed by Luterbacher *et al.* [28], we stress that these winter atmospheric data fields are assumed to be correct over the period of the reconstruction. Subsequently, we apply the same technique to obtain the Hilbert transformed versions of the proxy indices. However, these contain unresolved long time-scale trends longer than the training period, which hamper our effort to reconstruct the wave climate based on the six-term model. To obtain the long time scale in our reconstruction, we make the following assumption.

By assuming the long time-scale variation is sufficiently captured by three (proxy) climate indices, we can perform the reconstruction using the three-term model as previously applied in [12]. Based on the same assumption, we can also perform the reconstruction using a modified six-term model, containing the short time-scale fluctuations (less than 50 years) from the six-term model with the long time-scale variation (more than 50 years) described by the three-term model. Figure 6*a* shows the reconstructed one in 100 year *H*_{mp} at St Kilda, Haltenbanken and Valhall from 1661 to 2012, obtained by a combination of proxy indices from 1661 to 1961 and known indices from 1962 to 2012, using the modified six-term predictor model (solid lines) and the three-term model (dashed lines). Overall, the reconstructed results are similar for each location, and during the training period, the agreement between the two models is very close. Also shown on the same figure is the normalized histogram of the reconstruction for each location, these are approximately normally distributed. When the reconstructed results are plotted against each other over the period of reconstruction, the similarity between the two models becomes more evident, demonstrating the robustness of the two models in capturing the short time-scale fluctuations, as shown in figure 6*b*. However, because the six-term model is able to explain the short-term variability better than the three-term, the reconstruction from the modified six-term model should be more reliable.

It is worth remarking that Luterbacher's earlier data were generally based on terrestrial observations and this represents the marine region less well [29]. It was only from approximately 1850s onward that marine data well into the Atlantic were included. Thus, there has to be uncertainty in the reconstructed results in particular the long-term trend prior to the 1850s, and hence our reconstructions back then may not be as reliable as more recently. From figure 6*a*, we observe that in the 1990s the extreme wave height reaches the highest level ever seen (compared with the historic record) which aligns with the strongly positive phase of the NAO, and in the period from 1960 to 1990, there is a faster increase rate of the extremes for all locations. Also, there appear to be long-term trends at Haltenbanken and Valhall that may be statistically significant but these are influenced by apparently low activity prior to 1850s and by the recent heightened activity during the early 1990s.

Using the same reconstruction analysis, figure 7*a* shows the reconstructed one in 100 year *H*_{sMax} for all three locations, using the modified six-term predictor model. It can be seen that the fluctuations at all time scales are very similar to those of the reconstructed one in 100 year *H*_{mp}, which becomes more evident when the reconstructed results of *H*_{sMax} and *H*_{mp} are plotted against each other over the period of reconstruction. The solid red line is a 2:1 line showing the relative magnitude of *H*_{mp}∼2×*H*_{sMax}. This is in contrast to the usual estimate of *H*_{mp}=1.86×*H*_{sMax} assuming a 3 h sea-state with *N*∼1000 waves. The factor of ∼2 implies storms last longer than 3 h, perhaps 3× as long on average, with *N*∼3000 waves. This is one indication that *H*_{mp} is a more robust measure of storm severity than *H*_{sMax} [14,15,17].

Figure 7*b* shows the reconstructed 99th percentile of *H*_{s} for the same locations using the modified six-term predictor model. It can be seen that the variability of the reconstructed 99th percentile of *H*_{s} is similar to that of the 100 year *H*_{mp}, reflected by the *R*^{2} measured relative to the 100 year *H*_{mp}. The close similarity between the two reconstructions demonstrates the robustness of our extreme value analysis, given that the 99th percentile requires no extrapolation and hence it is more robust measure of extremes (though it is not directly informative in a return period or reliability framework for offshore structural design for example). It should also be noted that these values of the 99th percentile of *H*_{s} are substantially smaller than any of the threshold *H*_{mp} values used for the POT analysis, hence we are observing not just teleconnections influencing the extremes, but also more values further down into the bulk of the distributions.

It is also of interest to reconstruct mean wave climate, represented by annual mean *H*_{s} for comparison with the reconstructed extreme wave climate. For such climates, we apply similar correlation methodology described in [12] using the three-term predictor model (with additional high-pass filtering applied to the EA and the SCA indices), and the six-term predictor model as introduced in this paper (without any filtering). The reconstruction of the mean *H*_{s} is then performed using the modified six-term predictor model based on the same assumption, and a 5 year moving average is run across the reconstructed signal. The reconstructed mean *H*_{s} for the same locations is shown in the same figure 7*c*. It is worth noting that because of the filtering applied in the three-term model for mean *H*_{s}, the length of the reconstructed mean *H*_{s} (from 1668 to 2005) is slightly shorter than that of the one in 100 year *H*_{mp}. The long-term mean and standard deviation of the reconstructions for all locations are provided in table 3.

From the wave climate reconstruction, similar temporal variation is observed between the reconstructed extreme, 99th percentile and mean wave climate for most of the open North Atlantic locations (locations 1–4), over the same 5 year sliding window. In contrast, the temporal structures between the mean and extreme wave climates are different for the North Sea and the rest of the locations. The similarity is measured in terms of *R*^{2} relative to the 100 year *H*_{mp}, shown in the same figure and summarized in table 3 for all locations. Thus, the temporal variations in the extreme, 99th percentile and mean wave climates in these open North Atlantic locations are comparable, all these measures of the wave climates are influenced by the NAO because of the close proximity to the dominant storm tracks. For the rest of the locations including the North Sea, the extreme and mean wave climates are rather different: the mean wave climate is correlated to the NAO, but there is less evidence for the extreme wave climate. Also on the same figure, it is worth noting that the variation in the 99th percentile of *H*_{s} is closer to that of the 100 year *H*_{mp} compared with the mean *H*_{s}, demonstrating that the use of the 99th percentile as an in-sample estimate of extremes is a reasonable approach.

Overall, our reconstructions show an unusual increasing trend in both extreme and mean wave climates during the period from the 1960s to 1990s. This is consistent with [31,30], who reported an increasing trend in mean wave climate in this region, and [2], who conducted extreme value analysis on *H*_{sMax} for locations close to Valhall, Bruce and Haltenbanken. However, our long time-scale analysis implies that the increasing trend during that period was strongly influenced by the variability of the ocean–atmosphere system driven mostly by the NAO. Our reconstruction seems to suggest that the wave model and wave observations of both *H*_{sMax} and *H*_{mp} during the strongly positive phase of the NAO in the early 1990s should be satisfactory for determining design wave specifications for coastal and offshore structures going forward, though the observations from around 1900 would be equally applicable. In general, large variations in the extreme wave climate are observed, which should be accounted for in any long-term survivability analysis for coastal, marine and offshore activities over time scales of several decades. Unfortunately, analysis of the autocorrelation structure suggests that there is little internal predictability, as shown in figure 8. A rapid drop of the predictability is observed within 5 years. Nevertheless, whether the next winter is rough or less rough might be important for offshore operations (a drilling campaign for example) as opposed to platform survivability.

## 6. Conclusion

Long-term extrapolation of extreme wave height based on hindcast wave data is presented, using a POT technique for locations in the northeast Atlantic and the North Sea. The estimates of the extreme wave height using a simple exponential fit to the entire 54 years of record are consistent with the storm severity each location experienced. We present our estimates based on the most probable maximum individual wave height in a storm or *H*_{mp}, which is a measure of storm severity. We compare the results of the one in 100 year *H*_{mp} with the one in 100 year *H*_{sMax} as well as the in-sample 99th percentile of *H*_{s}, and conclude that our estimate of *H*_{mp} is robust. To account for possible non-stationarity in extreme value analysis, we perform 5 year sliding window analysis to produce variability of the estimated extreme wave heights in time, which we then correlate with the NAO and other atmospheric modes. Improved correlation is obtained using the six-term predictor model compared with the three-term model during the period of available data, with the six-term model including both the climate indices and their Hilbert transforms. The long-term variability less than 50 years (the length of the wave hindcast) is correlated to the decadal structures of the NAO and other modes, whereas some of the intradecadal fluctuation over less than 5 years is probably owing to finite size sample variability in the POT analysis. Overall, we find that the westerly extreme wave heights at all locations both in the open North Atlantic and in the North Sea are correlated to the NAO, and stronger correlation is observed for locations in close proximity to the dominant storm tracks.

With the two predictor models, we produce two types of reconstructions of extreme wave climate at each location, assuming the large-scale trend longer than the period of available data is sufficiently captured by the NAO and other indices. Both reconstructions are robust in describing the short timescale fluctuations (less than 50 years), however, because the six-term model is able to explain the variability slightly better than the three-term model, the reconstruction from the modified six-term should be more reliable. Our reconstruction shows in general large decadal variability of extreme wave height for all locations. In particular, there is a recent heightened activity of extremes in the 1990s, which aligns with the strongly positive phase of the NAO, and there is a faster increase rate of the extremes for all locations in the period from 1960 to 1990, though there has to be uncertainty in our reconstructed results particularly the long-term trend prior to 1850s where Luterbacher's pressure reconstruction may not be as reliable as more recently. Reconstruction of mean wave climate in terms of annual mean *H*_{s} is also performed using similar methodology. Comparison with the historic extreme wave climate reveals that the variability in time of the mean, 99th percentile and extreme wave climates for locations in close proximity to the dominant storm tracks in the open North Atlantic are comparable, whereas the wave climates for the rest of the locations including the North Sea are rather different. From the reconstruction of 100 year return period, we observe that *H*_{mp}∼2×*H*_{sMax}, which implies that storms typically last longer than 3 h, perhaps 3× as long on average.

In this paper, we have concentrated on the 100 year wave (formally the 10^{−2} yr^{−1} probability level) and its decadal variability. This could be a realistic design criterion for a wave power installation, with a design life of a few decades and the consequences of failure being just financial. The design criteria for manned offshore platforms will be much more severe, with (much) longer return periods being used. In general, the historic extreme wave analysis reveals high levels of decadal variability for all locations, both in the open North Atlantic and also in the North Sea. Analysis of the autocorrelation structure suggests little predictability, unless forecasting of the NAO into the future can be improved. The recent results of Smith *et al.* [32] are encouraging in this regard.

## Ethics

We declare we have no ethical concerns.

## Data accessibility

The climate indices are available from the NOAA/CPC at ftp://ftp.cpc.ncep.noaa.gov/wd52dg/data/indices/tele_index.nh. The reconstructed 500 mbar pressure maps by Luterbacher *et al*. [28] are available at https://www.ncdc.noaa.gov/paleo/pubs/luterbacher2002/luterbacher2002.html. Please contact the corresponding author for the proxy indices.

## Authors' contributions

H.S. and P.H.T. developed the methods used in this paper, which were implemented by H.S. All authors helped analyse the data and contributed to writing the paper.

## Competing interests

We declare we have no competing interests.

## Funding

This work is supported by UK Engineering and Physical Sciences Research Council (EPSRC; project EP/J010316/1 Supergen MARine TechnologY challenge), whose support is gratefully acknowledged.

## Acknowledgements

We thank BP Sunbury for providing the wave data, and Dr Tim Woollings of AOPP in Oxford for key inputs to this work and helpful discussions. We are also grateful for the help of two Oxford undergraduates Victoria Barker and David Bishop for preliminary work.

- Received May 20, 2016.
- Accepted August 17, 2016.

- © 2016 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.