Return levels of extreme European windstorms, their dependency on the North Atlantic Oscillation, and potential future risks

Windstorms are the most damaging natural hazard across western Europe. Risk modellers are limited by the observational data record to only ∼ 60 years of comprehensive reanalysis data that are dominated by considerable inter-annual variability. This makes estimating return periods of rare events difficult and sensitive to the choice of the historical period used. This study proposes a novel statistical method for estimating wind gusts across Europe based on observed windstorm footprints. A good description of extreme wind speeds is obtained by assuming that gust speed peaks over threshold are distributed exponentially, i.e. a generalised Pareto distribution having a zero shape parameter. The threshold and tail scale parameter are estimated at each location and used to calculate estimates of the 10- and 200-year return levels. The North Atlantic Oscillation (NAO) is particularly important for modulating lower return levels and modulating the threshold, with a less detectable influence on rarer extremes and the tail scale parameter. The length of historical data required to have the lowest error in estimating return levels is quantified using both observed and simulated time series of the historical NAO. For reducing errors in estimating 200-year return levels of an independent 10-year period, a data catalogue of at least 20 years is required. For lower return levels the NAO has a stronger influence on estimated return levels, and so there is more variability in estimates. Using theoretical estimates of future NAO states, return levels are largely outside the historical uncertainty, indicating significant increases in risk potential from windstorms in the next 100 years. Our method presents a framework for assessing high-return-period events across a range of hazards without the additional complexities of a full catastrophe model.

1 Introduction

Extratropical cyclones (ETCs) are the dominant weather system across Europe in the winter season, contributing most of the precipitation (Hawcroft et al., 2012) and strong winds (Ulbrich et al., 2001). ETCs are commonly associated with extreme wind speeds (also known as windstorms) that can impact infrastructure, agriculture, and transport and cause loss of life (e.g. Browning, 2004; Schwierz et al., 2010; Kendon and McCarthy, 2015). The strongest storms such as Lothar (26 December 1999) and Kyrill (18 January 2007) caused total insured losses exceeding EUR 5.5 billion (2022 equivalent; PERILS AG, 2022), with economic losses far exceeding this.

Although long-range predictability of European windstorms exists (e.g. Scaife et al., 2014; Befort et al., 2019; Degenhardt et al., 2023), extreme events are hard to predict and are simulated based on historical events and known natural variability. The main tools used to estimate rare events and their risk potential are catastrophe models (Grossi and Kunreuther, 2005). These models generally use reanalysis or climate model data as the driving hazard along with varying exposure and vulnerability data to quantify these rare risks. Catastrophe models must be able to quantify losses at long return periods, for example at the 200-year return level in order to comply with the Solvency II directive¹. Complying with this solvency directive becomes more difficult as there are only a few decades of coherent observational data from the historical period. Therefore, it is paramount to be able to understand and accurately quantify the risk of an event that would be at (or exceed) the 200-year return level. These risk estimates are often not freely available as catastrophe models are often licensed products with undisclosed vulnerability curves and calibration methods. Furthermore, catastrophe model event sets are often derived from century-long reanalyses or ensembles of climate model data, both of which may require substantial calibration or feature spurious trends (e.g. Bloomfield et al., 2018). Therefore, there is a clear need for a simple, transparent method to quantify high return levels associated with European windstorms for the wider risk community.

Projected extremes are dependent on the underlying distribution, and it is therefore important to understand the total variability in extreme events in a historical period (Woo, 2019). Research on extreme precipitation has shown an increase in the intensity of potential extremes when considering a larger sample sizes (Thompson et al., 2017). This is particularly important for European windstorms due to the pronounced decadal variability in storm numbers (Donat et al., 2011b; Dawkins et al., 2016; Cusack, 2023) and losses (Klawa and Ulbrich, 2003). With a varying length of historical record this may result in particularly stormy periods, such as the early 1990s, having a stronger influence on projected risks. Shorter records that do not cover this period would likely underestimate the risk potential due to the absence of this key component of historical variability. Consequently, a key question that arises is how variable estimations of extreme gusts are with historical records of differing lengths.

A well-known driver of the frequency and strength of cyclones and windstorms in Europe is the North Atlantic Oscillation (NAO; Hurrell et al., 2003). The NAO drives these variations through modulations of the large-scale pressure gradient across the North Atlantic Ocean, with more positive NAO phases leading to more cyclones with higher intensities than neutral or negative phases (Pinto et al., 2009; Gómara et al., 2014; Dawkins et al., 2016). With a changing climate, it is expected that positive NAO states will become more common (Fabiano et al., 2021) and cyclone wind speeds will become stronger (Priestley and Catto, 2022). It should therefore be expected that the NAO would modulate extreme gusts across Europe in both a historical and future climate, although the extent of this has currently not been explored.

Consequently, the four main questions to be addressed in this study are as follows:

1. Can 200-year return level gust speeds from European windstorms be reliably estimated using openly available observed windstorm footprints?

2. How important is the NAO in modulating high-return-level gust speeds?

3. How does the length of the historical record and choice of period influence return level estimates?

4. How could expected changes in the NAO under climate change lead to changing return levels?

2 Data

2.1 Footprint database

For analysing historical windstorms, a set of spatially coherent and validated footprints are required at high spatial resolution. Windstorm footprints are summaries of the overall hazard and are maps of maximum wind gust speed over a period of 72 h. The Windstorm Information Service (WISC; WISC, 2017) project provided a set of footprints based on historical events for analysing the range, severity, and impact of historical windstorms (e.g. Koks and Haer., 2020; Welker et al., 2021). Footprints were produced by the UK Met Office through the downscaling of ERA-Interim and ERA-20C reanalysis data. The ERA-20C footprints range from 1940–2010 and ERA-Interim footprints from 1979–2014. In this study, the ERA-20C footprints from 1950–1979 and the ERA-Interim footprints from 1980–2014 are used. There are some differences in footprints generated from the two reanalyses, although robust differences in footprint intensity cannot be determined. The footprints provide spatial estimates at 4.4 km resolution of the maximum 3 s wind gust at 10 m elevation over the 72 h period of each downscaled storm (WISC, 2017). The WISC footprints are similar to those provided by the XWS project (Roberts et al., 2014) but extend further back in time and are downscaled to a higher spatial resolution.

Figure 1 Storm footprints from the WISC catalogue. (a) Footprint of Storm Daria on 26 January 1990. Average footprint of two different length time periods: (b) 1972–2013 and (c) 1989–2013. (d) The difference between the two periods as a percentage. Units are metres per second (m s⁻¹).

The storms to be downscaled were selected from a set of objectively identified cyclone tracks (Hodges, 1994, 1995), with those chosen being of relevance to the insurance industry and also those exceeding a pre-determined intensity threshold (Steptoe, 2017). The time period of the cyclone track for downscaling is defined as ± 36 h of the maximum intensity of the storm, with the maximum intensity defined as the time of maximum 925 hPa wind speed over land. The downscaling is performed in two steps. Firstly the reanalysis data acts as boundary conditions for a 12 km version of the Met Office Unified Model that is initialised on 4 consecutive days for 30 h. Following this, the 12 km model acts as boundary conditions for the 4.4 km model. The 4 consecutive days of 4.4 km downscaled data is then concatenated and the relevant time period of the track extracted. The final footprint is then the maximum wind gust value at each grid point. A spatial Gaussian filter is applied to the footprint in order to eliminate spurious extremes generated in the downscaling. All the footprints have been validated against observations and result in 124 high-resolution windstorm footprints for the period 1950–2014 for use in this analysis.

2.2 NAO data

The NAO data used in this study are obtained from the NOAA CPC (https://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/nao.shtml, last access: 12 December 2023) and are calculated utilising rotated principal component analysis as described in Barnston and Livezey (1987). NAO indices are calculated using rotated principal component analysis applied to monthly 500 mb geopotential height fields from 1950–2000 that are then interpolated to daily values. NAO daily values are available for the period 1950–present and are standardised using the mean and standard deviation of the 1950–2000 climatological values. The NAO index for each storm is taken from the middle day of the 72 h WISC footprints.

Figure 2 Tests of the Gumbel fit for high return periods at the three representative locations: (a) Bergen (60.4^∘ N, 5.1^∘ E), (b) London (51.5^∘ N, 0.4^∘ W), and (c) Madrid (40.4^∘ N, 3.8^∘ W). The gusts for each location from 1950–2014 are plotted against their return period. The dotted black line indicates the 0.7 quantile, the red line indicates the Gumbel fit, the dark yellow line is the 10-year return period, and the blue line is the 200-year return period. The 95 % confidence intervals are estimated from a gamma distribution with shape and scale defined from the σ parameter and number of footprints and are shown by the grey shaded region.

2.3 Statistical methodology

For building the statistical model to estimate high-return-period windstorm gusts at any given location or grid point across Europe, the return period of gusts from the WISC footprints (e.g. storm Daria; Fig. 1a) first needs to be quantified. The return period is the expected time between successive events and therefore is equal to the reciprocal of the rate of the process. Therefore, the formulation of the return period (in years) of a wind speed (y) for an event Y>y is

$$\begin{array}{}\text{(1)}& T\left(y\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }\mathrm{P}r(Y>y)}},\end{array}$$

where λ is the rate of storm occurrence (WISC footprints per year; ∼ 2.1). Windstorms can be separated into events with strong winds (S) and those without strong winds ($\stackrel{\mathrm{‾}}{S}$). These events have different rates, λ_S and ${\mathit{\lambda }}_{\stackrel{\mathrm{‾}}{S}}$. For large return periods, it is assumed that the rate of occurrence for strong windstorms is considerably larger than for weak windstorms (${\mathit{\lambda }}_{\stackrel{\mathrm{‾}}{S}}\mathrm{Pr}(Y>y\mid \stackrel{\mathrm{‾}}{S})\ll {\mathit{\lambda }}_S\mathrm{Pr}(Y>y\mid S)$). As the WISC data were constructed just for meteorologically intense and impactful windstorms, we define this as our strong storm subset, and therefore we know the rate of strong storms with strong gusts, thus following these assumptions $T\left(y\right)\approx \left({\mathit{\lambda }}_S\mathrm{P}r\right(Y>y\mid S){)}^{-\mathrm{1}}$.

The exceedance probability can then be factored as follows:

$$\begin{array}{}\text{(2)}& \mathrm{Pr}(Y>y\mid S)=\mathrm{6}\mathrm{Pr}(Y>y\mid (Y>u)\cap S)\times \mathrm{Pr}(Y>u\mid S),\end{array}$$

where u is a threshold that is large enough to use extreme value theory generalised Pareto fits to the model $\mathrm{Pr}(Y>y\mid (Y>u)\cap S)$. By assuming that the tail shape parameter is zero and in the Gumbel domain of attraction, this gives exponentially distributed excesses above the threshold:

$$\begin{array}{}\text{(3)}& \mathrm{Pr}(Y>y\mid (Y>u)\cap S)=e^{\left(-\frac{y-u}{\mathit{\sigma }}\right)}.\end{array}$$

In Eq. (3), σ is the tail scale parameter that can be estimated by method of moments (Bowman and Shenton, 2006) from the mean excess above the threshold (u), yielding

$$\begin{array}{}\text{(4)}& \widehat{\mathit{\sigma }}={\displaystyle \frac{{\sum }_{y_i>u}{y}_i-u}{{\sum }_{y_i>u}\mathrm{1}}}.\end{array}$$

Furthermore, the quantity $p\left(u\right)=\mathrm{P}r(Y>u\mid S)$ can be estimated from the relative frequency of exceedances:

$$\begin{array}{}\text{(5)}& \widehat{p}\left(u\right)={\displaystyle \frac{\mathrm{1}}{n}}\sum _{y_i>u}\mathrm{1},\end{array}$$

and so if u is taken to be the qth empirical quantile (u=y_q), then $\widehat{p}\left(u\right)=\mathrm{1}-q$. As the WISC data are already a sub-selection of extreme wind speeds, a relatively low quantile (q=0.7) is chosen to model the tail of the distribution. The sensitivity of the estimated return levels to the value of u are discussed below.

Equating $\mathrm{Pr}(Y>y|S)=({\mathit{\lambda }}_{S}T{)}^{-\mathrm{1}}$ from Eq. (1) with $\mathrm{Pr}(Y>y|S)=\exp (-(y-u)/\mathit{\sigma }\left)p\right(u)$ from Eqs. (2) and (3) gives a simple prediction for the T-year return level that is based on the qth empirical quantile (u), the mean excess above threshold ($\widehat{\mathit{\sigma }}$), and the rate of windstorm occurrence (${\widehat{\mathit{\lambda }}}_S$):

$$\begin{array}{}\text{(6)}& \widehat{y}=u+\widehat{\mathit{\sigma }}\left(\text{log}T+\text{log}\widehat{p}\left(u\right)+\text{log}{\widehat{\mathit{\lambda }}}_S\right).\end{array}$$

In order to test the functionality of the method, the statistical model is applied to the WISC data at three example locations. These three locations vary in latitude and represent varying influence from the leading pattern of variability, the NAO. The three locations are Bergen (60.4^∘ N, 5.1^∘ E), London (51.5^∘ N, 0.4^∘ W), and Madrid (40.4^∘ N, 3.8^∘ W). Figure 2 shows a demonstration of the statistical model. The three locations have different gust speed distributions, with higher gusts at the northernmost locations. The estimation of the 10-year return level at each location is very similar from our model and from the empirical data. However, the 200-year return levels are larger than any gust from the WISC footprints at each location. For Bergen the 200-year estimate is 45.9 m s⁻¹ (39.7 m s⁻¹ WISC max), for London 39.0 m s⁻¹ (38.6 m s⁻¹), and for Madrid 30.1 m s⁻¹ (24.1 m s⁻¹). The statistical model has also been applied to the five largest cities across Europe (not shown) and to a set of windstorm footprints derived from the ERA5 re-analysis (Appendix A) and provided similar results. Our model is therefore able to estimate high-return-period gusts that are unprecedented in the observations.

Figure 3 Hazard maps generated by the statistical model for the (a) 10-year and (b) 200-year return levels across Europe. (c) The 0.7 quantile threshold used to calculate the return levels. (d) The mean excess (sigma) used in the return level calculations. All units are metres per second (m s⁻¹).

The chosen empirical quantile is 0.7 for estimating the extreme gusts. Figure S1a–c in the Supplement demonstrates how the 200-year return level varies with a changing value of u (or q). When u is very small, there are large variations in (and large values of) the 200-year return level due to the inclusion of low gust speeds in our model fitting. There is a stabilisation in the 200-year return level for q≥ 0.7, with variations within 5 m s⁻¹ for all three locations. This is the case up until the very highest thresholds, where larger variations in the 200-year return level are again noted; however, these high thresholds have a reduced number of footprints in the model fitting, so as with the low thresholds, these are unlikely to provide realistic estimates of the high return levels.

3 Results

3.1 Wind gust hazard maps

By applying this method to all grid points of the WISC footprints, estimates of the 200-year return levels can be produced for all of Europe (Fig. 3). The 10-year return level is also examined as it is a return level that is easier to validate with historical claims data. The 10-year (Fig. 3a) and 200-year (Fig. 3b) return levels show a maximum in the return levels across northwestern Europe and the North Atlantic Ocean. This is expected as this is the region of strongest gust speeds in the WISC footprints (Fig. 1b and c). Return levels decrease radially from the UK, with lower values across Iberia, Italy, eastern Scandinavia, and eastern and southeastern Europe. These lower values are associated with a lower frequency and reduced intensity of windstorms in the historical catalogue due to a greater distance from the main North Atlantic storm track (Priestley et al., 2020). The largest 200-year return levels exceed 60 m s⁻¹ to the northwest of Ireland. Over land, the largest 200-year return levels are 55–60 m s⁻¹ across eastern Scotland, northeastern England, and Denmark, with values above 40 m s⁻¹ for the majority of the rest of the British Isles, northern France, Benelux, and Germany (Fig. 3b).

Figure 4 Maps of the (a) ${\widehat{\mathit{\beta }}}_{\mathrm{0}}$, (b) ${\widehat{\mathit{\beta }}}_{\mathrm{1}}$, (d) ${\widehat{\mathit{\alpha }}}_{\mathrm{0}}$, and (e) ${\widehat{\mathit{\alpha }}}_{\mathrm{1}}$ parameters for regressions applied to the WISC footprints from 1950–2014 with daily NAO covariate. Panels (c) and (f) show p values of the $\widehat{\mathit{\beta }}$ and $\widehat{\mathit{\alpha }}$ regressions respectively. Units are (a, d) metres per second (m s⁻¹) and (b, e) metres per second per SD ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{SD}}^{-\mathrm{1}}$).

The dominant parameters in the return level estimates are the model threshold (u; Fig. 3c) and the mean excess above the threshold ($\widehat{\mathit{\sigma }}$; Fig. 3d). The spatial pattern of u (Fig. 3c) is similar to the return level maps (Fig. 3a and b) in that it features a maximum in the northwest of the domain in the area surrounding the North Atlantic. The magnitude of u decreases toward southern and eastern Europe. This pattern suggests a strong NAO influence, with larger values in the region where the NAO drives high wind speeds and damaging windstorms (Pinto et al., 2009). Unlike u, $\widehat{\mathit{\sigma }}$ does not have the same north–south dipole over western Europe (Fig. 3d). Instead values vary much less, with only a 4–5 m s⁻¹ variation across the entire WISC domain. The variation in $\widehat{\mathit{\sigma }}$ is not coherent and is therefore consistent with this component being independent of large-scale patterns, such as the NAO.

3.2 Generalising the model to include an NAO covariate

The parameter u displays a pattern that is indicative of an NAO influence and can therefore be generalised to include such variations. If the threshold is chosen to be the qth empirical quantile, then $u={\widehat{y}}_q$ and then ${\widehat{y}}_q$ and σ can be generalised to be a function of a covariate x. The covariate x will be the daily NAO index at the time of a WISC footprint occurrence. By estimating the threshold using linear quantile regression and estimating σ using generalised linear regression for a gamma distribution with identity link, this yields

$$\begin{array}{}\text{(7)}& {\displaystyle }{\widehat{y}}_q={\widehat{\mathit{\beta }}}_{\mathrm{0}}+{\widehat{\mathit{\beta }}}_{\mathrm{1}}x,\text{(8)}& {\displaystyle }\mathit{\sigma }={\widehat{\mathit{\alpha }}}_{\mathrm{0}}+{\widehat{\mathit{\alpha }}}_{\mathrm{1}}x.\end{array}$$

Figure 5 Dependence of gusts on the NAO at the three representative grid points: (a) Bergen (60.4^∘ N, 5.1^∘ E), (b) London (51.5^∘ N, 0.4^∘ W), and (c) Madrid (40.4^∘ N, 3.8^∘ W). Scatter points are the WISC grid point gusts against the daily NAO index. Horizontal dashed black line indicates the 0.7 quantile for each location. The red line is the q=0.7 estimation based on the NAO gust regression. The dark yellow line indicates the 10-year return level, and the blue line is the 200-year return level using the Gumbel domain estimation.

Applying these regressions to the WISC data yields the parameters shown in Fig. 4. The $\widehat{\mathit{\beta }}$ parameters (Fig. 4a and b) exhibit the typical NAO-influence pattern, with larger values and positive regression coefficients across the northwest of the WISC domain and smaller values and negative coefficients across eastern, southeastern, and southern Europe. This indicates that more northerly latitudes and northwestern European locations have a positive NAO–gust threshold relationship and vice versa. Furthermore, p values of this regression (Fig. 4c) indicate widespread significance of the NAO influence for northwestern and northern Europe, with large areas of p<0.001. The $\widehat{\mathit{\alpha }}$ parameters (Fig. 4d and e) do not feature an NAO influence, and values are more variable across all of Europe. The north–south dipole is not evident, and it therefore appears that the NAO has minimal influence on the windstorm gust excesses. The significance testing (Fig. 4f) indicates there are isolated areas where the NAO influences σ, yet these are not widespread or indicative of a mode of variability. Despite this finding for the $\widehat{\mathit{\alpha }}$ parameters, it may be that an NAO relationship across larger areas of Europe is undetectable in our small data sample and that with a larger pool of windstorm footprints a signal may emerge. Due to the undetectable signal in the $\widehat{\mathit{\alpha }}$ parameters the NAO covariate is only applied to the threshold (u) in the statistical model, and therefore Eq. (6) becomes

$$\begin{array}{}\text{(9)}& \widehat{y}={\widehat{\mathit{\beta }}}_{\mathrm{0}}+{\widehat{\mathit{\beta }}}_{\mathrm{1}}x+\widehat{\mathit{\sigma }}(\text{log}T+\text{log}(\mathrm{1}-q)+\text{log}{\widehat{\mathit{\lambda }}}_S).\end{array}$$

Figure 5 shows the NAO covariate model (Eq. 9) tested at the same three grid points as previously. Each location features a different NAO relationship. The most northerly point, Bergen, has a strong positive relationship between return level and the NAO (p<0.01). An NAO of +1.5 standard deviations yields a 200-year return level of > 50 m s⁻¹. For the London grid point, this NAO relationship is still positive, although weaker than for Bergen (p< 0.1), as would be expected. Finally, for the Madrid grid point an insignificant (p> 0.5) negative–neutral relationship is present, and therefore when the NAO is more positive, the estimated return level gets smaller. The reduced NAO sensitivity for the Madrid grid point is to be expected from the in-land and southerly nature of this location.

Figure 6 NAO dependence of return levels. Maps of the (a, c) 10-year and (b, d) 200-year return level across Europe using the NAO covariate. Maps are shown for an NAO signal of (a, b) +0.5 and (c, d) −0.5. Units are metres per second (m s⁻¹).

Applying the NAO covariate method to all grid points yields the return levels for the 10- and 200-year return periods shown in Fig. 6. The same regional variation and distribution in return levels across Europe are seen as in Fig. 3. Return levels are higher for the 200-year return period (Fig. 6b) than the 10-year (Fig. 6a) with values ∼ 10 m s⁻¹ higher across northwestern Europe. It should be noted that the calculated 200-year return levels are higher when using the NAO covariate (Figs. 6b and 3b), with maximum values of > 65 m s⁻¹ over the North Sea and North Atlantic Ocean and values exceeding 50 m s⁻¹ across the northern British Isles, Denmark, and northern Germany. The variability in return level with an NAO state of +0.5 (Fig. 6a and b) and −0.5 (Fig. 6c and d) is also apparent for both return levels, with higher values of up to 5 m s⁻¹ across the northwest of the domain for the NAO+ state (see also Fig. S2).

The spatial pattern of the 10- and 200-year return levels (Fig. 6) differs due to the differing contribution of u and the excesses ($\widehat{\mathit{\sigma }}\left(\text{log}T\right)$) with return period (Fig. S3). The contribution of the NAO-dependent u decreases with increasing return period, and the contribution of the NAO-independent $\widehat{\mathit{\sigma }}\left(\text{log}T\right)$ increases. Therefore, as the influence of the NAO is only detectable in the regression of u, it can be said that the relative importance of the NAO on our estimated return levels decreases with an increasing return period.

3.3 Sensitivity to choice of historical period

One factor that can contribute to the uncertainty in the estimation of the 200-year return level is the length of the historical catalogue. This is especially important for (re-)insurers and their need to understand risk in the next 10 years, as this is a time horizon for business planning, and 10 years is also the typical maximum service life of a catastrophe model. The varying length of catalogue has a substantial impact on the average footprint (Fig. 1b–d), and when applied to the statistical model at a return period of 200 years, these differences are likely to be amplified. Therefore, the understanding of future risk is likely to differ with these different length catalogues. The 200-year return level estimates vary considerably when using historical catalogues of different lengths (Fig. S4). The mean squared error (MSE) of the 200-year return level for historical catalogues from a year in the range of 1951–2014 relative to the full catalogue of 1950–2014 increases with shorter catalogue lengths (Fig. S4a) with grid point return level estimates being both under- or overestimated depending on the historical catalogue used (Fig. S4b and c).

In order to quantify how different catalogue lengths compare in their return level estimates, an independent time period is required for validation. For this, a 10-year time period of WISC footprints is taken (e.g. January 2005–December 2014), and the 200-year return levels are estimated. The historical catalogues for the preceding years starting at a length of 1 year and extending all the way out to 55 years (i.e. 2004, 2003–2004, 2002–2004, …, 1950–2004 in this example case) are taken and return levels estimated. The MSE of all catalogue return levels are then calculated relative to the 10-year validation period to quantify the optimal historical catalogue length required to minimise the error in estimated return levels. To account for natural variability all possible 10-year validation periods are tested with start years ranging from 1955 to 2005.

Figure 7 Mean square error in the 200-year return level estimation of historical catalogues of different lengths against a subsequent 10-year period from the WISC catalogue for (a) Bergen, (b) London, and (c) Madrid. The solid black line shows the median mean squared error when using all possible periods. The dark and light grey areas represent the 50 % and 95 % confidence interval on the standard error respectively from the catalogues with different validation periods.

Applying this methodology to the same three grid points noted above we find that all locations (Fig. 7) feature a high MSE for the shortest catalogue lengths. The high MSE with short catalogue length is expected due to the limited sample of footprints contributing to the return level estimate. Increasing catalogue lengths to above 2 years causes a sharp reduction in the MSE at all locations, with limited improvements after ∼ 10–15 years. Therefore, historical catalogues longer than 10–15 years do not yield improvements in the return level MSE at these three locations. There are variations in the MSE estimate from year to year and large uncertainty in the MSE, which is an artefact of the limited data catalogue. With a limited historical catalogue, potential reductions in uncertainty or MSE for catalogue lengths longer than 40 years cannot be quantified. Greater uncertainty in MSE is also seen for the longest catalogue lengths at Madrid (Fig. 7c), which is a sampling artefact from having fewer long catalogue samples from the WISC data. Validation periods of different lengths (i.e. 5, 20 years) were also tested with conclusions being insensitive to this change (not shown).

3.4 Simulating the NAO for a more robust catalogue length estimation

For a more confident estimation of the optimal catalogue length and to assess the impact of data catalogues longer than 40 years in length, a time series of the NAO and the windstorm model parameters are simulated. The NAO is simulated as a sinusoidal variation of a period of 70 years (e.g. Wanner et al., 2001; Olsen et al., 2012; Wang et al., 2017) and amplitude of 1.5 standard deviations to match the variations from the WISC footprints (not shown). The WISC footprints occur at a rate of ∼ 2.1 per year, and the model is built to fit the top 30 % of these events; therefore, occurrences of extreme storms and their NAO phase can be estimated from these frequencies. This entire simulation is done for 1000 years (Fig. 8a).

Using the regressional relationships between the WISC footprints and the NAO state (Fig. 4a and b), the value of u (Fig. 8a) and σ can be estimated in the simulated time series and the resultant gust determined (Fig. 8b). Comparing the WISC grid point gusts for Bergen with our simulation (Figs. 5a and 8c), it is notable that the simulated gusts exceed those in the WISC dataset considerably, with gusts in excess of 60 m s⁻¹. This is due to the longer time series in our simulation compared to the WISC dataset.

Figure 8 Summary of simulated events at Bergen based on 1000-year NAO time series. (a) NAO phase and simulated storm threshold (u) at time of simulated events. (b) Simulated event gusts (m s⁻¹). (c) NAO phase plotted against the simulated gusts. The solid red line shows the q=0.7 regression line.

Figure 9 Mean square error in the (a) 10-year and (b) 200-year return level estimation of historical catalogues of different lengths against a subsequent 10-year period with events from a 1000-year simulation for Bergen. The solid black line shows the median MSE from all possible periods. The dark and light grey areas represent the 50 % and 95 % confidence interval on the standard error respectively from the catalogues with different validation periods. Vertical dashed grey lines indicate the periodicity of the NAO used in simulations.

Using this simulated time series, the 10- and 200-year return levels can be estimated for different catalogue lengths for our three locations, with analysis performed as in Fig. 7. At Bergen the highest MSE is seen for the shortest catalogue lengths (Fig. 9). As in Fig. 7a, this is a result of the large inter-annual variability and the chance that 1 year will be representative of the following 10 being very small. For the 200-year return levels (Fig. 9b), a large reduction in the MSE is then seen as catalogue length increases to 15 years. For catalogue lengths of 15–50 years there are minimal changes in the MSE. Following this there is a slight reduction in MSE associated with capturing the full cycle of the imposed NAO cycle; however, there is overlap in the confidence intervals at 20 years and 70 years. Beyond 70 years and unlike the WISC results, there is now little variation in the MSE out to catalogue lengths in excess of 250 years and a much reduced spread and variation in the median. Therefore, there is limited benefit to using historical catalogues longer than ∼ 20 years when estimating the 200-year return level; however, if possible, a catalogue spanning a full NAO cycle is preferable. It should be noted that in Fig. 9 the MSE never reaches zero, indicating that a perfect return level estimate is never achieved. This is a result of unresolved natural variability that cannot be captured regardless of catalogue length.

At the 10-year return level (Fig. 9a) a reduction in MSE occurs up to a catalogue length of ∼ 10 years. Unlike the 200-year MSE, there is an oscillation in the MSE that has a period of ∼ 70 years and results in a maximum MSE for catalogue lengths of 40 years and a minimum after ∼ 80 years. The ∼ 70-year oscillation is consistent with the defined periodicity of the NAO simulation. The oscillation is present as the NAO influence is more pronounced for shorter return periods (Fig. S3). Consequently, the 10-year return levels more strongly reflect the variation in NAO, whereas the 200-year return levels will not. To reduce the MSE, the NAO signal of the historical catalogue relative to the 10-year validation period is important. For shorter catalogue lengths the NAO signal (and simulated gusts) will be similar; however, catalogue lengths of ∼ 35 years feature an NAO signal that is most different from the 10-year validation period and hence the MSE will be most different. When a full NAO cycle is sampled (∼ 70 years), the full proportion of NAO variability is captured and there is a minimum in MSE. For catalogue lengths longer than the prescribed NAO periodicity, the average NAO signal varies less and a dampening of the MSE variation is seen.

The structure of MSE evolution seen at Bergen (Fig. 9) is similar at London (Fig. S5) and Madrid (Fig. S6). For the 10-year return levels the oscillatory evolution at Bergen (Fig. 9a) is less apparent due to the weaker influence of the NAO on return levels (Figs. 4 and 5). As a result it is less important to sample a full cycle of the NAO, and the MSE varies little for catalogue lengths longer than ∼ 20 years.

3.5 Change in return levels from theoretical future NAO states

Recent estimates of the historical trend of the NAO have been upwards at a rate of +0.15 standard deviations per decade (1950–2020; Blackport and Fyfe, 2022). If it is assumed that the historically averaged NAO state is neutral, then a theoretical state of +1.5 standard deviations in 100 years is a plausible upper estimate. For estimating return levels of gusts in this theoretical future NAO state, events are simulated. The role of natural variability on return levels is also considered, and to account for this a 70-year period is used that features two 35-year periods at $+(-)\mathrm{0.5}$ standard deviations compared to the reference state. For the historical state the reference is NAO neutral, and events are simulated at −0.5 and +0.5. For the future state the reference NAO is +1.5, and events are simulated at +1 and +2. The 70-year periods are simulated 100 000 times for our three grid points to obtain robust return levels. Selecting a model threshold and estimating the return levels of this two-component, variable NAO state model differs slightly from the earlier methodology and is documented in Appendix B.

Figure 10 Boxplots of estimated (a–c) 10-year and (d–f) 200-year return levels for (a, d) Bergen, (b, e) London, and (c, f) Madrid from simulated gusts. Return levels are estimated from 100 000 70-year simulations with NAO phase that varies to ± 0.5σ of a set NAO state. Boxplots show the median return level of these 100 000 simulations, and boxes extend to the 25th and 75th percentiles, with whiskers extending to the 2.5th and 97.5th percentiles. Red crosses signify the theoretical return level from a two-state NAO. Dashed horizontal grey lines indicate the 95 % confidence interval of the WISC 10- and 200-year return levels.

Using simulated gusts from our variable state model, the estimated return levels at Bergen for the NAO state of +1.5 are largely outside the WISC uncertainty for the 200-year return level (Fig. 10d). Therefore, if historical trends were to continue for the next 100 years, then the 200-year return levels would likely be unprecedented in this theoretical state. This is more evident at the 10-year return level (Fig. 10a) due to the increased influence of the NAO on these shorter return periods. For all NAO states the return levels estimated using the simulated, variable events are positively biased relative to theoretical values using exact parameters (Fig. 10). This is a result of the additional natural variability present in the simulated events from the varying NAO state (see Appendix B). Therefore, models that assume a stationary historical climate are also likely to be positively biased relative to actual return levels.

Similar changes in return levels are also seen for the future NAO state at London (Fig. 10b and e). However, as the NAO has a reduced influence the increase in 200-year return levels is less apparent with most estimates within the historical uncertainty (Fig. 10e). Nevertheless, the potential increases in 10-year return levels are still large enough to be greater than the historical uncertainty (Fig. 10b). At Madrid (Fig. 10c and f), the NAO influence is so small that any changes in the return levels are insignificant compared to the historical uncertainty.

4 Conclusions

This study has presented a method for estimating high-return-period gusts for the European domain based upon a set of observed windstorm footprints. This simple and transparent method is agnostic to the input dataset and operates without requiring calibration or the complexities of a full catastrophe model. However, it does not replace catastrophe models as no loss estimation or vulnerability module is included. The amount of historical data required to minimise errors in the return level estimation has been tested, and the sensitivity of return levels to theoretical future NAO states has been quantified. We raised several questions in Sect. 1, and the key conclusions that answer these are as follows:

• Return levels at the 200-year return period have been quantified and associated uncertainties quantified. Return levels are higher over northwestern Europe, with the NAO playing a role in increasing (decreasing) return levels across northern (southern) Europe.

• The NAO is important for setting the threshold for gust speed return levels, with the excess above the threshold being stochastically generated. The NAO modulates the lower return levels and is key in determining the threshold, with particular significance over northwestern Europe. The tail scale parameter is much more variable with a less detectable influence from the NAO.

• Only 20 years of historical data is required to reduce errors in 200-year return level estimates made from the independent and subsequent 10-year period. Increasing historical records longer than this do not lead to meaningful improvements in return level estimates due to unresolved natural/inter-annual variability. Estimates of 10-year return levels are more sensitive to the length of historical catalogue due to the larger NAO influence.

• Under theoretical future NAO states the 10- and 200-year return levels are likely to be unprecedented relative to historical values. Return level estimates are outside the historical uncertainty across northern Europe, with minimal changes across southern Europe. This would lead to considerably higher gusts and impacts than those currently considered.

Our modelling approach makes several simplifying assumptions. Firstly, for estimating the tail a zero shape parameter is used. This implies that an infinite return level is possible for the longest return periods. This is inaccurate, and in reality the shape parameter would be negative as there are frictional processes, instabilities, and limits to the amount of transferrable kinetic energy from a windstorm. However, we are confident in our estimations as the generated return curves are a good fit for the input data for numerous locations across Europe. These results provide additional data points to complement estimations made from more complex catastrophe models that will help guide risk modellers and aid understanding of potential European windstorm risks.

With regards to the data used, there is a historical limit due to the constraints of observations of severe European storms. This limit to 1950 may result in considerably large events that would constrain the model being missed or key components of the natural variability being excluded. However, with ongoing efforts to accurately represent historical storms (Hawkins et al., 2022), this may be improved in the future. Alongside this, the data used from the WISC project are only a small sample of all windstorms over Europe representing just the most intense events, and this may lead to a skewing of our return level estimations toward the worst case scenario.

Our methodology allows for the assessment of theoretical future return levels without the need of climate change simulations. However, these theoretical return level estimations only represent the changes in the circulation aspects of the NAO. One factor not considered in the large increase in return level for northern latitudes of Europe is the thermodynamic contribution. As a large amount of the potential increase in cyclone wind speeds across Europe is a result of moist processes (Dolores-Tesillos et al., 2022; Binder et al., 2023), it is likely that gusts would increase further beyond those predicted in the model here. These changes could have further implications on potential impacts. Only the NAO was considered a low-frequency modulator of windstorm intensity in our model. However, other phenomena have been shown to influence European winter weather through connections to the NAO. These are phenomena such as tropical precipitation (Scaife et al., 2017, 2019), the Atlantic multi-decadal oscillation (Börgel et al., 2020), and the El Niño–Southern Oscillation (Zhang et al., 2019a, b) to name a few. These low-frequency phenomena all have partial connection to the NAO and therefore may play a role in modulating windstorm strength through either the threshold or excesses in this model. Investigating potential linkages further may lead to improvements in the model and therefore more accurate return level estimations.

There are numerous future applications of our model. Firstly, this model could be used to investigate future changes in return levels. With projected future changes in the NAO and other leading modes of European weather variability (Fabiano et al., 2021; Blackport and Fyfe, 2022), the simulation of the NAO can be modified to apply trends and changes in variability to determine how this would affect the higher-return-period gusts. This could provide vital insight to the projected end-of-century increase in European storminess (Donat et al., 2011a; Zappa et al., 2013; Priestley and Catto, 2022; Dolores-Tesillos et al., 2022). Furthermore, seasonal forecasts of the NAO could be imposed to inform decision makers at shorter timescales. Consideration should be taken, however, due to the unresolved natural variability that is present in return level estimates regardless of the length of historical record. Finally, despite this model being developed for European windstorm risks, it could be easily adapted for other hazards/risks and act as a framework for assessing extreme events of other perils. Validation and testing would be required to meet regulatory and business requirements, but this model can act as an open-source and transparent method for investigating and quantifying a range of natural hazards.

Appendix A: Applying the statistical model to ERA5 footprints

In testing the applicability of the statistical model it has also been applied to a set of footprints from the ERA5 reanalysis (Hersbach et al., 2020). Footprints were created from a set of extratropical cyclone tracks. Tracks were identified following the method of Hodges (1994, 1995), which uses 850 hPa relative vorticity as a tracking variable. Tracks were created from 1-hourly data for the entire calendar year (1 January–31 December) for the period of 1980–2020. Full details of the cyclone identification and tracking method can be found in Hoskins and Hodges (2002). Tracks were filtered to those passing through a defined North Atlantic–European domain (30–75^∘ N, 40^∘ W–40^∘ E) to encapsulate all tracks that would pass through continental Europe and to cover all the WISC domain.

Figure A1 Wind gust footprint from ERA5 of track forming on 20 January 1990 and dissipating on 28 January 1990. Coloured shading represents the maximum grid point wind gust in metres per second (m s⁻¹). White line and crosses indicates the track and hourly track points respectively.

Footprints were created using hourly near-surface (10 m) wind gusts. These data provide the maximum 3 s wind gust at each grid point in the previous hour. At each track time step the gusts in the surrounding 5^∘ are associated with the storm, and the maximum gust at each grid point across the lifecycle of the track is retained in order to create a coherent footprint. This method closely follows that of Roberts et al. (2014). An example footprint is shown in Fig. A1.

These footprints are then applied to the statistical model described in Sect. 2.3, and the analysis is performed as in Fig. 2. As there are many more footprints in our ERA5 dataset than from WISC, there are a considerably greater number of points to fit our model (Fig. A2a–c). As a result of the increase in the number of footprints, which are mainly of lower gust speeds, a higher threshold quantile for ERA5 is used than for the WISC footprints, with q=0.9 for the ERA5 footprints. Despite the increase in the number of footprints, the higher value of q, and the coarser resolution of the ERA5 data, the model fit is a good approximation for the data, and estimates of the 200-year return levels are able to be made.

The estimated 200-year return levels (Fig. A2a–c) are different from those estimated from WISC (Fig. 2a–c), as would be expected due to the differing input data. However, the values are consistent with those from WISC, with the main difference being an increase in the value of the 200-year return level for the London grid point (Fig. A2b). What is further re-assuring about the applicability of our model is the stability of the 10- and 200-year return levels for high quantiles used in the model (Fig. A2d–f). Estimations of these two return levels vary by less than 5 m s⁻¹ for q≥0.9.

Figure A2 (a–c) As Fig. 2 but using ERA5 footprints as an input. (d–f) The sensitivity of the 10-year (solid black) and 200-year (solid red) return levels to the choice of quantile threshold for (d) Bergen, (e) London, and (f) Madrid. Vertical black lines indicate the 0.5 and 0.9 quantiles. Thin red and black lines represent the 95 % confidence interval based upon the estimate of the σ parameter.

Appendix B: The role of natural variability

To gain insight into the effect of natural variations in modulators of the extremes, it is useful to consider a simple two-component mixture model. Consider a hazard process that is an equal mixture of two processes Y_a(t) and Y_b(t) that have exponential exceedances above thresholds a and b>a respectively. Then for y≥b,

$$\begin{array}{rl}\mathrm{Pr}(Y>y)& =\mathrm{0.5}\mathrm{Pr}(Y_a>y|Y>a\left)\mathrm{Pr}\right(Y_a>a)\\ & +\mathrm{0.5}\mathrm{Pr}(Y_b>y|Y>b\left)\mathrm{Pr}\right(Y_b>b)\\ & =\mathrm{0.5}\text{exp}(-(y-a)/\mathit{\sigma })\mathrm{Pr}(Y_a>a)\\ & +\mathrm{0.5}\text{exp}(-(y-b)/\mathit{\sigma })\mathrm{Pr}(Y_b>b).\end{array}$$

If thresholds a and b are chosen to be the qth quantiles of Y_a and Y_b, then

$$\begin{array}{rl}\mathrm{Pr}(Y>y)& =\mathrm{0.5}\text{exp}(-(y-a)/\mathit{\sigma })(\mathrm{1}-q)\\ & +\mathrm{0.5}\text{exp}(-(y-b)/\mathit{\sigma })(\mathrm{1}-q)\\ & =\mathrm{0.5}(\mathrm{1}-q)({\text{e}}^{a/\mathit{\sigma }}+{\text{e}}^{b/\mathit{\sigma }}){\text{e}}^{-y/\mathit{\sigma }}.\end{array}$$

Hence, the T-year return level is given by

$$y_T=\mathit{\sigma }\text{log}\left(\mathrm{0.5}\right({\text{e}}^{a/\mathit{\sigma }}+{\text{e}}^{b/\mathit{\sigma }}\left)\right)+\mathit{\sigma }(\text{log}T+\mathrm{1}-q+\text{log}{\mathit{\lambda }}_S),$$

which equals that derived in Eq. (7) if one sets $u=\mathit{\sigma }\text{log}\left(\mathrm{0.5}\right({\text{e}}^{a/\mathit{\sigma }}+{\text{e}}^{b/\mathit{\sigma }}\left)\right)\in [a,b]$. By integration of the probability density function, it is straightforward to derive the mean excess of this process above threshold $u\in [a,b]$:

$$E(Y-u|Y>u)=\mathit{\sigma }+{\displaystyle \frac{b-u}{{\text{e}}^{-(u-a)/\mathit{\sigma }}+\mathrm{1}}},$$

which can be seen to be greater than or equal to σ. Hence, estimating σ by assuming it is equal to the mean excess will lead to a positive bias that vanishes only as u→b. In summary, the effect of natural variability in the tail location parameter is to cause the mean excess to overestimate the tail scale parameter, which in turn will lead to overestimates of return values for longer return periods. Put simply, if natural variability is not taken into account, the magnitudes of the extreme events will appear to be more variable than they actually are.