Evaluating the Effects of Disjunct Sampling and Averaging Time on Extreme Wind Speeds

1. Introduction

Extreme wind speeds for the design of structures are typically assessed by extreme-value analysis (EVA) of wind speeds observed at 10m above ground level in flat open country and embodied into codes of practice and other design guidance, e.g. as in [1]. The accuracy of these empirical assessments depends on many factors, including the characteristics of the anemometers, their exposure in terms of upwind roughness changes and their density over the geographical region. More recently, synoptic-scale wind speeds have been deduced at the surface from the atmospheric pressure fields of reanalysis models [2], again at hourly intervals, developed from the pioneering study by [3] – a methodology that resolves issues of geographical cover.

Most anemometers are situated at airports, reporting mean and maximum gust speed and direction over a 10-minute period at hourly or, latterly, ½-hourly intervals to the World Meteorological Organization (WMO) standards. Reports are archived for more than 35,000 stations around the globe in the US National Centers for Environmental Information (NCEI) Global Hourly, Integrated Surface Data (ISD), dataset DS3505 1. These standard reports are under-sampled (disjoint), as they omit the wind speeds occurring in the 50, or 20, minutes between the reports. Some meteorological authorities also make contiguous observations, e.g. mean and maximum gust values over a one-hour period at hourly intervals, so that no gust speeds are lost but it is still possible to miss a higher average wind speed that spans across intervals. The WMO standards (FM12 SYNOP and FM15 METAR) report wind speed in integer knots (1kn = 0.5144ms^-1) and directions in 10° increments but, to accommodate speeds in other units, DS3505 archives speeds in integer 0.1ms^-1 increments (e.g. 20kn is archived as 103). Contiguous hourly-mean observations for the UK are available in the Met Office MIDAS Open database hosted by the Centre for Environmental Data Analysis 2. Starting in 1998 at a few stations and implemented at more than 800 stations by 2004, the US Automated Surface Observation System (ASOS) provides observations of the 2' mean and the maximum 5" gust speed and direction at 1' intervals in the ISD dataset DS6405 3 and is currently the only readily accessible dataset sampled at such a high rate.

National and international codes of practice differ in their definition of extreme wind speeds for the design of structures. The UK was the first to move to a probabilistic basis with the introduction of a 3-second gust speed with a 50-year return period [4] on the basis that gust speeds were thought to be less affected by ground roughness than mean speeds. However, extreme gust speeds are much more likely than mean speeds to include high-value artefacts and just one artefact per year can dominate the EVA [5]. With the development of an atmospheric boundary layer model for heterogeneous terrain [6], allowing reliable exposure corrections for mean speeds, the UK moved to hourly-mean speeds applied with a gust factor model for gust speeds [7]. Development of the Eurocodes led to the adoption of 10-minute mean wind speeds [8] as the datum which aligns with the WMO observational standards. Each of these different definitions remain in force somewhere around the globe.

This present study focuses on two principal issues: a) the effect of averaging period and b) the effect of sampling interval on the assessment of extreme mean wind speeds, specifically on $\overline{V}\{\overline{t},R\}$, the $R$-year return speed averaged over time $\overline{t}$. Resolving the first issue enables transformation of values between averaging periods and resolving the second issue allows underestimates of peak speeds to be corrected. The essence of these issues is illustrated in Figure 1 which presents typical 5h traces of: (a) $\overline{t}=10\text{'}$ and (b) $\overline{t}=60\text{'}$, running mean wind speed at $t_{∆}=1\text{'}$ intervals by the ◯ (circle) symbols. Note that wind speed units are maintained in the original observed WMO standard of knots (1kn = 0.5144ms^-1) throughout this study to avoid unit conversion bias [5]. For the $\overline{t}=10\text{'}$ running mean in (a), contiguous sampling at $t_{∆}=10\text{'}$ intervals is shown by the thin (red) line, and disjoint sampling at $t_{∆}=60\text{'}$ intervals by the thick (blue) line. It is evident that the disjoint samples miss the central peak speed as well as the short-duration (mesoscale) peak at minute 63, while the contiguous samples slightly underestimate these peaks. For the $\overline{t}=60\text{'}$ running mean in (b), the short-duration peak is suppressed by the running mean and contiguous sampling at $t_{∆}=60\text{'}$ intervals underestimates the central peak by a larger amount. In sampling a continuous timeseries the probability of missing a peak is controlled by the relative timescale, defined here as the disjoint sampling ratio:

$$∆=t_{∆}/\overline{t}=\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{a}\mathrm{l}/\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{d}$$

(1)

Hence $∆>1$ represents disjoint sampling, $∆=1$ contiguous sampling, and $∆<1$ over-sampling. For the ASOS DS6405 $1\text{'}$-interval wind speeds, $∆(\overline{t}=10\text{'})=1/10$ and $∆(\overline{t}=60\mathrm{\text{'}})=1/60$, which are the closest to continuous that can be achieved with the ASOS observations.

The remainder of this paper is organized as follows: §2 is a review and discussion of relevant previous studies; informing §3 which describes the models, methods and data used in the two stages of this study; §4 presents the results of each stage with relevant discussion, and finally; §5 assembles and discusses the conclusions of the study which are summarized by the bullet points in §6.

2. Review of Previous Studies

Previous studies of these issues are reviewed chronologically to take account of insights and flaws that are addressed, or not, by later studies:

1. In [9] a new wind map was proposed for Italy for $\overline{V}\{R=50\mathrm{y},\overline{t}=10\text{'}\}$ to the Eurocode standard. This was based on EVA of $\overline{V}\left\{\overline{t}={10}^{\text{'}}\right\}$ for 42 stations at $t_{∆}={180}^{\text{'}}$ intervals and for 27 stations at $t_{∆}=10\text{'}$ intervals, corrected for exposure. Preliminary studies suggested that the 3-hour interval underestimates in the range 5% to 10%;

2. In [1], proposing a new zonal wind map for Germany, a factor of $f_{\overline{t}}=1.06$ was applied to transform from hourly-mean observations, $\overline{t}=60\text{'}$, into the Eurocode standard ten-minute mean, $\overline{t}=10\text{'}$, without revealing the source of this value. The study included thunderstorm downbursts and gust fronts in addition to the synoptic scale windstorms, fitting only the (least reliable) far-right tail of the EVA distribution of this mixture to evaluate the dominant component in any given location. The study proposed an optimum sample size, $N_{opt}$, defining the rank $m=N_{opt}$ of lowest observation as given by $P\left\{m=N_{opt}\right\}=1-P\left\{m=1\right\}$. This resulted in fewer fitted observations than recommended by [10] for the XIMIS EVA method, but was necessary to isolate the dominant component of the mixed wind climate;

3. The [3] study was the first to derive extreme wind speeds at the surface from maps of sea-level atmospheric pressure and was the inspiration for the rapid development of this technique into using re-analysis models, e.g., [2]. 489 individual storms with a central pressure below 970mbar were selected to represent the most intense storms affecting Europe between 1953 and 1995. Geostrophic wind speeds in $60’$ time steps over a 0.5° interval grid covering Europe were evaluated by linear interpolation between the 6h-interval maps over the lifetime of each storm. 50-year return geostrophic speeds were estimated by the peaks-over-threshold (POT) method and transformed to surface wind speeds using the geostrophic drag law. The principal advantage of this approach is that the surface wind speed is predicted everywhere in the field to the same resolution, which is particularly valuable where anemometer stations are sparse;

4. Study [11] focused on the principal issues of the present study: empirically through observations and theoretically when modelled as a Gaussian Markov chain. Empirical assessments of averaged speeds at 5 locations in Denmark and 3 locations in the Arabian Gulf were made in terms of the mean of the observed annual maxima. Attenuation factors were evaluated for sample intervals from $t_{∆}={10}^{\text{'}}$ to $t_{∆}=12\mathrm{h}$ and for averaging periods, $\overline{t}$, over the same range – apparently assuming the two effects to be mutually independent. Evaluating their empirical fits to the relevant averaging times and sample intervals indicates that the factor for $\overline{t}={10}^{\text{'}}$ from $\overline{t}={60}^{\text{'}}$ is $f_{\overline{t}}=1.054$, compared with 1.06 from [1], and that $t_{∆}={60}^{\text{'}}$ sample intervals underestimate yearly maxima by 3.1% and $t_{∆}={180}^{\text{'}}$ intervals by 6.0%. The principal deficiency of the Markov model is that each new value depends only on the previous value, whereas wind speeds have a memory of previous values over several hours, as expressed by the integral timescale, $T$, of their autocorrelation$;$

5. Study [12] addressed three issues: 1 – The starting wind speed of anemometers; 2 – The appropriate extreme-value distribution; and 3 – The averaging time and sampling interval between observations. The first issue was demonstrated to have no significant effect on extremes and, in any case, the starting speed of the latest generation of fast-response cup-anemometers is very much improved and is effectively zero for sonic anemometers. While debate over the second issue continues, the argument against asymptotic models and in favor of the sub-asymptotic, or “penultimate” form of the Fisher-Tippet Type 1 (FT1), or Gumbel distribution is demonstrated by the latter’s effectiveness in the XIMIS method [10,13]. Hence the only third issue of averaging time and sampling interval is relevant to the present context. The empirical analysis was similar to [11], starting with $\overline{t}=10\text{'}$ mean wind speeds at hourly intervals, re-sampling at longer intervals, $t_{∆}=$2, 3, 4, 6, 8, 12 and 24h, then extrapolating a regression fit back to contiguous $t_{∆}=10\text{'}$ intervals. They claimed that $t_{∆}=10\text{'}$ intervals represent the theoretical limit value, but this neglected the underestimation of continuous running averages by contiguous averages illustrated in Figure 1. They concluded that the correction factor to contiguous $t_{∆}=10\text{'}$ average speeds from $t_{∆}=$ ${30}^{\text{'}}$, ${60}^{\text{'}}$ and ${180}^{\text{'}}$ intervals was $f_{∆}=$ 1.03, 1.05 and 1.08, respectively;

6. Study [14] was the first to exploit the $t_{∆}=1\text{'}$ interval observations in DS6405, comparing hourly means from $t_{∆}=1\text{'}$ interval, $\overline{t}=2\text{'}$ average, DS6405 data to the mean from DS3505 at hourly intervals at 8 US stations. It claimed an 8% – 38% overestimate of $\overline{V}\{R=50\mathrm{y},\overline{t}=60\text{'}\}$ by DS3505 “due to large discrepancies” in DS3505. The DS3505 observations for these stations are typical hourly reports compiled from the DS6405 data over the $10\text{'}$ averaging period, between ${20}^{\text{'}}$ and $10\text{'}$ before each hour, so are ${\overline{t}=10}^{\text{'}}$ averages and not $\overline{t}=60\text{'}$ averages, as was assumed. One principal contribution to the overestimates is therefore the missing $\overline{t}={60}^{\text{'}}{10}^{\text{'}}$ transformation factor $f_{\overline{t}}\approx 1.06$ [1]. Another is that the ${\overline{t}=10}^{\text{'}}$ average includes mesoscale convective events which are suppressed by the $\overline{t}=60\text{'}$ average. The short observation period of 14 years leads to inherent overestimation of $\overline{V}\{R=50\mathrm{y},\overline{t}\}$ by EVA due to statistical variance. These effects do not explain the largest overestimates, but EVA is notoriously susceptible to any high-valued artefacts in the record. The “extensive QA/QC process” of this study was not explained but needed to have been rigorous, as in [5]. Additional procedures have since been developed to identify and correct frequent artefacts particular to the ASOS DS6405 data [15]. The study gave the first estimate for the slight increase in the rolling hourly average (continuous) over the block hourly average (contiguous) of 2% – 3%;

7. Study [16] used Monte-Carlo simulations of Weibull parents by the [17] method, using the probability density and spectrum obtained from DS6405 data at each of 6 stations in central USA to obtain simulated records of 18,000y because “the combined effect of the sample size and disjunct sampling cannot be well understood by only using one set of data with limited years of observation”. The locations were chosen to avoid hurricane contributions, and the observations were stripped of thunderstorm events indicated by the DS3505 “present weather” codes. The study claimed that the Gaussian Markov chain model in [11] overpredicted the effect of disjunct sampling for ${\overline{t}=10}^{\text{'}}$ averages at large sample intervals. It predicted that the factor from $t_{∆}=60\text{'}$ intervals to contiguous $t_{∆}=10\text{'}$ intervals is $f_{∆}=1.04$, and to $t_{∆}=1\text{'}$ intervals is $f_{∆}=1.06$ (where the latter can be regarded to be effectively a continuous running average). The study concluded that the effect of disjunct sampling depends only on the disjoint ratio, $∆$ (Equation 1);

8. Study [18] reviewed the previous studies, above, and concluded that none were “satisfactory and reliable”, being particularly critical of the Gaussian process and zero-memory Markov assumptions in [11] and of the accuracy of the extrapolation to contiguous sampling in [12]. The study used contiguous ${\overline{t}=10}^{\text{'}}$ average wind speeds at 117 stations across Germany to calibrate a closed-form equation relating yearly maxima, $\overline{V}\{R=1\mathrm{y},\overline{t}=10\text{'}\}$, to disjoint samples, confined specifically to synoptic storms, of the form:

$$f_{∆}={\left[1-a\times \mathrm{l}\mathrm{n}(∆)\right]}^{-1/b}\mathrm{f}\mathrm{o}\mathrm{r}∆\ge 1$$

(2)

where $f_{∆}$ is the correction factor, and $a,b$ are empirical values to be fitted. Two issues with Equation (2) are that the datum for $f_{∆}=1$ was $∆=1$, contiguous not continuous, and the boundary conditions as $∆\to 0$ (continuous) are not met. In calibrating Equation (2) it was assumed that the ${\overline{t}=10}^{\text{'}}$ average eliminated contributions from thunderstorm events, but their influence is only reduced. The fitted values had no geographical correlation, leading to $a=0.156$ and $b=4.5$ as being applicable to the whole of Germany, with $\sigma \left(f_{∆}\right)=0.0114\times \log ∆$. Other useful values emerged from the data, including the Weibull shape factor $w=1.82\pm 0.24$ which is comparable to UK values [19];

The chronologically evolving estimates of the transformation factors $f_{\overline{t}}$ and $f_{∆}$ for $\overline{V}\{R=50\mathrm{y},\overline{t}=10\text{'}\}$, summarised in Table 1, show consistent trends, but inconsistent values. Having critically reviewing earlier work, the values of $f_{∆}$ proposed by [18] are considerably larger. A third relevant timescale is absent from all the previous studies: the integral timescale, $T$, derived from the autocorrelation, ${\rho }_V\left\{\delta t\right\}$, by integration or by fitting an appropriate model. This may account for the variation in values for the same $∆$ between locations in different wind climates.

Most of the previous studies focus on synoptic-scale events. By including the annual and daily harmonics in the spectra, the simulation in [16] by the [17] spectral method raises a fundamental issue. Spectra are formed by an integral-average through time which also removes all phase information, so that simulations generated using [17] method are forced to be stationary and so the harmonics are resolved homogeneously through time as cycles of constant amplitude and unknown phase. The phase of the fundamental harmonics can be deduced from the physics, but not the Oort and Taylor [20] sidebands that represent the seasonal modulation of the diurnal cycles. This creates anomalies in the simulated timeseries, including negative velocities at night in periods of relative calm, as the diurnal cycles have the same amplitude at all wind speeds. This is incompatible with the current understanding that each wind climate is caused by several meteorological mechanisms which, when examined in a seasonal-diurnal analysis, resolve into a disjoint mixture of phase-locked, amplitude-modulated random components [21,22]. By working in the time domain, so preserving phase, only the Harris methodology is currently capable of simulating the seasonal-diurnal components of a disjoint mixture [23].

3. Materials and Methods

The present study addressed the issues empirically by analyzing observations, and parametrically by analyzing long-period simulations, using the ASOS parameters: ${\overline{t}=2}^{\text{'}}$ averages and $t_{∆}=1\text{'}$ sampling intervals as the datum, while taking due account of insights gathered from the earlier studies. All analyses of this study were implemented in the R statistical computing language, Version 4.4.1 (https://cran.r-project.org/).

3.1. ASOS Observations

The ASOS observations from 2000-2003 to 2023 at $1\text{'}$ intervals for the 642 well-exposed stations across the contiguous USA (CONUS) in Figure 2, were downloaded from the ISD dataset DS6405. The stations were all at airports and are identified by their unique International Civil Aviation Organization (ICAO) 4-letter codes, e.g. KBOS for Boston, Logan International Airport. Following the approach of [24] and earlier studies, observations that might have come from hurricanes were removed over the three-day period centered on landfall, initially just at stations in the landfall States. Owing their rarity, hurricanes and tropical cyclones are assessed separately from other synoptic-scale storms, in terms of their characteristics as a class and their likelihood of occurrence derived from observations of their tracks [25]. This methodology has been recently extended to assess the action of severe North Atlantic depressions across Europe [26]. Secondly, all the convective gust events >20kn and non-meteorological artefacts from the [27] dataset of classified convective gust events were removed over the ${60}^{\text{'}}$ period centered on the peak speed, so as to retain only the synoptic class. The later EVA analysis exposed potential non-synoptic artefacts remained at only 74 of the 642 stations, all of which were examined manually, and confirmed artefacts were removed as reported in §3.4.

The standard ${\overline{t}=10}^{\text{'}}$ and ${60}^{\text{'}}$ averaging was applied to the $t_{∆}=1\text{'}$ interval observations by a running block-average, after which the observations were re-sampled at $t_{∆}={10}^{\text{'}},{30}^{\text{'}}$, ${60}^{\text{'}}$ and ${180}^{\text{'}}$ intervals. Re-sampling was repeated in subsets offset by $1\text{'}$ intervals to maximise the available information in the observations. For example, $t_{∆}=$ ${60}^{\text{'}}$ provided 60 subsets of observations for analysis, albeit mutually correlated.

3.2. Autocorrelations and Integral Timescales

Autocorrelations (non-dimensional) and autocovariances (dimensional) are tolerant of missing observations, whereas Fourier transformations for the spectra are not. Here the autocovariances were evaluated up to a maximum lag of 2 years after resampling at ${60}^{\text{'}}$ intervals to limit the computational overheads. Accounting for leap years by the average of 356.25 days per year, exactly 2 years ensured an integer number of ½-daily cycles and two whole annual cycles, ensuring that the integral timescales, $T$, were effectively the same with or without inclusion of the daily/annual harmonics.

Example results for KBOS, Boston, MA, are shown in Figure 3(a) for the first 10 days and (b) for the whole 2 years of lags. Figure 3(c) shows the autocovariance after removal of the harmonics by heterodyning [28], together with the best regression fits to the von Karman (vK) autocorrelation model, with its fixed 5/3 decay constant, and the bandwidth-limited pink noise (BLPN) model, introduced in [23] as a “universal replacement” for vK, since this allows for a variable decay constant and the limited bandwidth set by the averaging time, $\overline{t}$. Both models are good fits to the observed autocovariance and are virtually indistinguishable in the rapid decay towards the origin, justifying the use of the simpler vK model in §3.5. This provided three estimates for the integral timescale, $T$, which are mutually consistent, but around half the value found by [28] in the UK.

The integral timescale, $T$, is the scaling factor for time that normalizes the autocorrelation function with respect to time lag, and the corresponding spectrum with respect to frequency, into standard expressions. Neither of the other two time-dimensioned parameters, $O$ and $R$, are timescales in the conventional sense. The observation period, $O$, sets the population of independent peaks available for EVA, so controls the statistical variance of the results. The return period, $R$, is just a convenient proxy for the annual probability of exceedance, $\mathsf{\Phi }$, i.e., $\mathsf{\Phi }\left\{R=50\mathrm{y}\right\}=1-1/R=1-1/50=0.98$.

3.3. Weibull Analysis

The classic Weibull analysis of the observed wind speeds is relevant because the index, $w$, controls the rate of convergence of extremes to the FT1 distribution [29]. The scale $C$ and index $w$ parameters of the Weibull distribution:

$$P\left\{V\right\}=1-\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\left(V/C\right)}^w\right]$$

(3)

were estimated by regression after linearization to:

$$\ln \left(-\ln \left(1-P\right)\right)=w\times \ln \left(V\right)-w\times \mathrm{l}\mathrm{n}\left(C\right)$$

(4)

Following [30], $P\left\{V\right\}=P\left\{V<v\right|v>0\}$, is the cumulative probability of $V$ after removal of calms. The classical Weibull plot of $\ln \left(V\right)$ against $\ln \left(-\ln \left(1-P\right)\right)$ derives $w$ from the slope and $C$ from the intercept of a linear fit to Equation 4. The example in Figure 4(a) for Boston, MA, happens to be an excellent fit, helped by the 60-fold increase in observations at 1’ intervals as opposed to the more typical 60’ intervals, and the reduction in statistical variance this provides. But not all stations fitted quite this well. Figure 4(b) shows the probability density function (PDF) of $w$ for all stations at each of the averaging periods, with the corresponding ensemble mean, $〈w〉$, and standard deviation, $\sigma \left(w\right)$.

3.4. Extreme-Value Analysis

XIMIS [10] was adopted as the EVA method in this study because [13] recently demonstrated it to be the most reliable of the available methods for this context. XIMIS allows for incomplete convergence to the FT1 [29] distribution, for left-censoring of the peaks, for inverse variance weighting of the fit, and for preconditioning to maximize convergence to the penultimate FT1 model:

$$\mathsf{\Phi }\left\{y\right\}=\exp \left(-\exp \left(-y\right)\right)\mathrm{a}\mathrm{n}\mathrm{d}y=\left(V^w-U^w\right)/d^w$$

(5)

Here, $y$ is the reduced variate [29], $U$ is the mode, $d$ is the dispersion, and the convergence index $w$ is the index of the parent Weibull distribution [10].

Following experience over the last two decades in the use of XIMIS for supervised fitting, the left-censoring threshold was set to admit a maximum of 50 ranks (2½ times the observation period), but not more than half of the independent peaks. This allowed automated fitting to proceed unsupervised, which was essential given the huge number of fits required. The independent peaks in each subset were selected as the largest value occurring in the period between up/down-crossings of the mean speed, as illustrated in Figure 5, following the principle introduced by Davenport [31] in his adaption of the communication theory of Rice [32]. This predicts that there is one, and only one, statistically independent peak value between each successive up-crossing of the mean. As in [33], the running 100s mean was used to clearly define the mean-crossings. The Boston example in Figure 5 (a) illustrates the integer knot values of the reported $\overline{t}=2\text{'}$ mean, contrasting with the fractional values in (b) for the Rayleigh simulation. Figure 5(b) demonstrates the action of ensuring independence from adjacent peaks, showing that the secondary peaks either side of the highest are excluded because they occur between the same mean-crossings. The smaller peaks within their own mean-crossing periods on either side contribute, although they may subsequently be excluded by the left-censoring threshold of XIMIS.

XIMIS was applied at all ASOS stations for $\overline{t}=10\text{'}$ and ${60}^{\text{'}}$, using: a) the station values of $w$; b) $\overline{w}=1.7$ and 1.8, the means from Figure 4(b); and c) $w=2$, the Rayleigh value [34]. Figure 6 shows example XIMIS analyses for $\overline{t}=10\text{'}$ at KBOS, Boston, MA, as classical Gumbel plots, $V~y$. The ranked peaks for each fitted sub-set are shown by the + (red) symbols, the ensemble mean for each ranked subset by the ◯ symbols and their fit by XIMIS with $w=2$ by the thick curve. The outer chained curves are the 5%/95% confidence limits for the individual values, and the thin (red) curves are the 5%/95% confidence limits for values predicted from the fit. In the left-hand chart for $t_{∆}=1^{\text{'}}$ there is only one subset, so + and ◯ are coincident, but in the right-hand chart there are 60 subsets. The most volatile fits are for $\overline{t}=10\text{'}$ and $t_{∆}=1^{\text{'}}$because there is only one subset, so no ensemble to average. 74 stations (~12%) had outliers above the 95% confidence limit and, as recommended by Gumbel [29], these were individually checked. Most cases were found due to erroneous values, convective gusts <20kn, or events that evaded the quality control procedure, so that they were excluded from analysis. This left 20 valid single outliers (~3% of the stations) for $\overline{t}=10\text{'}$ due to short-duration (mesoscale) events, reducing them to 5 single outliers (~1% of the stations) for $\overline{t}={60}^{\text{'}}$as the additional averaging suppressed these events. These were retained to be addressed by the inverse variance weighting of XIMIS.

Ensemble averages were compiled of the mode and dispersion parameters in Equation (5) and of $\overline{V}\{R=50\mathrm{y}\}$ evaluated from the individual fits for each subset. As Equation (5) is linear in terms of each parameter value raised to the same index, $w$, it follows that $〈\overline{V}\{R,U,d\}〉\equiv \overline{V}\{R,〈U〉,〈d〉\}$, allowing $\overline{V}\left\{R\right\}$ to be evaluated for any other return period from $〈U〉$, $〈d〉$ and the associated $w$.

3.5. Simulation Method

The principal role of the second stage simulation was to examine the issues parametrically, using very long simulated time series to minimize statistical variance, and to evaluate how the parameters scale. Available simulation methods were the two time-domain methods of Harris: [35] for Rayleigh ($w=2$), [36] for Weibull with general index, and the more generalized spectral method of Nichols et al [17]. The choice need to keep the overall computation time within practical limits, requiring a compromise between best representation and speed. The Rayleigh method is the fastest, matching the target parameters in a single pass and its index, $w=2$, lies within the upper quartile for the ASOS stations, Figure 4(b). The computation time using either of the other methods would have been prohibitive. A secondary role of the simulation was to investigate the case of zero averaging period, $\overline{t}=0$, which is impossible from station observations as all include some degree of averaging.

Simulations of 1,000 years of wind speed were generated by the Rayleigh [35] simulation method at $t_{∆}=1\text{'}$ intervals for vK autocorrelations with a range of integral timescales, $T.$ This method exploits the fact that a Rayleigh timeseries, $w=2,$ may be formed as the orthogonal resultant of two Gaussian timeseries and that linear operations on Gaussian timeseries remain Gaussian, e.g., when applying an autoregressive (AR) filter. The required autocorrelation was evaluated, then transformed into AR filter coefficients by the Yule-Walker equations [37] which, when applied to the two Gaussian random timeseries, resulted in a Rayleigh-distributed timeseries with vK autocorrelation with the desired timescale, $T$. The same initial Gaussian streams were used to simulate each different timescale and averaging period to allow direct comparison between simulations and for compatibility with the empirical ASOS data. The simulation produced values normalized to unity mean.

Although the 1,000-year period was shorter than the 18,000 years used by [16], it proved quite adequate for the purpose, and allowed sub-division into independent observation periods, $O=10$y, 20y, 50y and 100y to emulate the range of observations found in practice, additional to the correlated sets. For example, with $O=20$y and $t_{∆}=60\text{'}$, typical of DS3505 data, there were 50 independent sets, each of 60 correlated subsets for EVA, requiring 3,000 XIMIS fits.

4. Results

4.1. ASOS Observations

4.1.1. Integral Timescales

The integral timescale, $T$, is the time-domain scaling factor for any given autocorrelation function. Figure 7 presents the relevant results of the ASOS autocorrelation analysis as color-mapped contour maps of CONUS, where: (a) for vK and (b) for BLPN show that the integral timescale, $T$, has a consistent geographical pattern with most variation in the high-topography western ⅓ of the USA and along the Gulf coast. Figure 3(c) showed an excellent match between vK and BLPN autocovariances for Boston, which is typical for all stations, with Figure 7(c) showing <10% difference overall and <2% difference over the eastern ⅔ of the USA. The BLPN decay constant, Figure 7(d) shows a similar pattern, but shows a gradient with latitude over the eastern ⅔.

Figure 8 investigates the BLPN decay index at stations east and west of longitude −100°. In (a) the PDFs are seen to separate into two distinct, but overlapping, distributions where the mean for the east stations is indistinguishable from the vK value of 5/3 but is significantly lower for the west stations. Note that each distribution shows a ripple matching the mode of the other, indicating some ‘leakage’ either side of longitude −100°. In (b), the quantile-quantile (Q-Q) plot shows a very small correlation between integral timescale and decay index, and that the difference between the all-station mean, 1.53, and the vK value is small compared with the range of scatter. Note that the four station values in the top left corner of (b) were taken to be invalid outliers and were excluded from the contour maps in Figure 7.

4.1.2. Weibull Parameters of Parent Wind Speeds

The Weibull index values are relevant to justify using the Rayleigh simulations with $w=2$. The PDFs in Figure 4(b) show that the ensemble mean, $〈w〉$ for all stations, increases with averaging period, $\overline{t}$. This is interpreted as progressive elimination of shorter duration events as $\overline{t}$ increases. The $\overline{t}=2\text{'}$ average includes contributions from the micrometeorological range of the atmospheric boundary layer (ABL) spectrum [38] which are progressively eliminated by the $\overline{t}=10\text{'}$ and $\overline{t}={60}^{\text{'}}$ averages which define the boundaries of the mesoscale range, the so-called Spectral Gap, so the residual mesoscale components for $\overline{t}=10\text{'}$ are eliminated by $\overline{t}={60}^{\text{'}}.$ Even so, with only synoptic-scale components, $〈w〉\{\overline{t}={60}^{\text{'}}\}=1.8$, which is still less than the ideal Rayleigh value predicted by Offset Elliptical Normal (OEN) model theory [19,39]. However, the synoptic-scale wind speed remains a mixture of disjoint components, typically cyclonic and anticyclonic systems. Figure 11 of [13] demonstrates how the index of a single Weibull distribution will present in the range $1.1\le w\le 2$ for a disjoint mixture of two Rayleigh components, depending on their relative frequency. The ranges of the PDFs in Figure 4(b) are much greater than the differences in $〈w〉$.

Figure 9(a) maps the geographical distribution of $w\{\overline{t}=10’\}$, showing lower values in the eastern and western ⅓ and higher values in the middle. The pattern for $w\{\overline{t}={60}^{\text{'}}\}$ is not shown because it is virtually identical (with 98% correlation and 0.042 standard error). In Figure 9(b), the Q-Q plot for $w\left\{\overline{t}=10’\right\}~w\{\overline{t}=60’\}$, the regression slope (solid line) matches the corresponding ratio of $〈w〉$ in Figure 4(b), and the dashed lines are the 5%/95% confidence limits on the standard error. The high correlation and small standard error indicate that this observed geographical variation is real.

4.1.3. Extreme Wind Speeds

The previous assessment of the reliability of EVA [13], showed XIMIS to be the most reliable when the index, $w$, in Equation (4) takes a fixed value, i.e. when using a two-parameter fit instead of a three-parameter fit, because fitting with $w$ as a free parameter is poorly conditioned. But the question remains, which value? The geographical variation of station values in Figure 9(a) is inconvenient when seeking a universal solution for mapping, and their mean values depend on averaging period, $\overline{t}$. The Rayleigh value $w=2$ is convenient in allowing direct comparisons with the Harris simulation model [35] of §3.2 but, more importantly, because that is predicted by the theory of the OEN model [19,39] and is expected to control the curvature of an OEN mixture distribution in the upper tail.

This question was resolved pragmatically by correlating $\overline{V}\{R=50\mathrm{y}\}$ for the individual station indices and the ensemble-mean for all stations, using the Rayleigh values as the datum in the Q-Q plots, Figure 10. In all cases the correlation coefficient, $\rho$, was unity to the fourth decimal place and the fitted slope is ~0.6% above unity. While the index, $w$, controls the curvature of the Gumbel plot, it also affects the mode and dispersion value which, for $O\approx 20\mathrm{y}$, rotates the fitted curve around a point close to $R=50\mathrm{y}$, making $\overline{V}\{R=50\mathrm{y}\}$ less sensitive to $w$.

Figure 11 maps the 1-minute interval, 10-minute and hourly average, 50y-return wind speeds across CONUS from the ASOS 1-minute interval observations, exposing the geographical variation. Although equivalent to design maps in codes and regulations, these maps are not intended to be used as such because the contributions of hurricanes are omitted and the “special zones” in the US code [40] are represented only when a station lies within the zone, e.g. KCTB, Cut Bank, MT for interior Chinook winds. However, Figure 11 bears no direct relationship to the design 3-second gust values in the US code [40] because these are dominated by thunderstorm downbursts over most of CONUS [41] and by gusts in hurricanes along the Gulf and eastern coasts [25] which are outside the scope of this study.

4.1.4. Disjoint Sampling Factor

The disjoint sampling factors, $f_{∆}$, of each of $t_{∆}={10}^{\text{'}},30\text{'},{60}^{\text{'}}$ and ${180}^{\text{'}}$ sampling intervals were estimated by using the ASOS $t_{∆}=1\text{'}$ sampling interval as the datum, i.e. $f_{∆}=\overline{V}\left\{t_{∆}\right\}/\overline{V}\{t_{∆}=1\}$, as proxy for true continuous data, leaving the consequences of this choice to be addressed later. Table 2 lists $〈f_{∆}〉$, the ensemble-mean values for all 642 ASOS stations, with their standard deviations, for the 50-year return period.

The expectation from [16] is that values of $〈f_{∆}〉$ for all averaging periods should collapse together in terms of the disjoint ratio, $∆$. Figure 12 tests this expectation for $\overline{t}=10’$ and $\overline{t}={60}^{\text{'}}$, in which the error bands indicate the 5%−95% confidence limits of the Gauss error function (erf) for each $〈f_{∆}〉$. Values of $〈f_{∆}〉$ decrease with increasing return period, while the error bands become wider. The collapse progresses from acceptable for $R=1\mathrm{y}$ in (a) to very poor for $R=1000\mathrm{y}$ in (c), indicating that $∆$ is not the sole controlling parameter and that other timescale ratios must be relevant, e.g. $\lambda =T/\overline{t}$. For comparison with [16], the thick (green) curve in each chart represents the model Equation (6) fitted to the ASOS data by optimization. This differs from Equation (2) in transforming to nominally continuous $t_{∆}=1’$, and converges properly to $f_{∆}=1$ as $∆\to 0$.

$$〈f_{∆}〉=1+A\times \mathrm{l}\mathrm{n}(1+B\times ∆)$$

(6)

Empirical values for the coefficients $A$ and $B$ are listed in Table 3. For the $t_{∆}$ intervals examined, Equation (6) gives $f_{∆}$ = 1.050 for ½-hourly, 1.072 for hourly, and 1.113 for 3-hourly observations. The [16] value from hourly to $1^{\text{'}}$ intervals is 1% lower, but is based on only 6 ASOS stations in a N-S line around −100° longitude.

Most of the earlier studies considered only the transformation to $t_{∆}=10\text{'}$ (contiguous), at which $f_{∆}$ = 1.026 for ½-hourly, 1.047 for hourly, and 1.088 for 3-hourly observations. The ${t_{∆}=180}^{\text{'}}$ value is in the upper quartile of the range suggested by [9]. The [11] values 1.03 for hourly and 1.06 for 3-hourly are both lower. The corresponding [12] values, 1.03, 1.05, and 1.08 respectively, match to the second decimal place. The [16] value for hourly is only 0.7% lower. The curve from the [18] formula for $f_{∆}$ cuts diagonally across this study’s model in Figure 12(b) close to the hourly interval, so their hourly value is 3% higher and the 3-hourly value 5% higher.

The model values in Figure 12 for $t_{∆}=1’$, at $∆=1/10$ for $\overline{t}=10’$ and at $∆=1/60$ for $\overline{t}={60}^{\text{'}}$ are first order estimates of the error of using $t_{∆}=1’$ as a proxy datum for true continuous data. These factors, in Table 4, could be used to adjust the empirical $〈f_{∆}〉$ estimates and to refit the model iteratively, but are too small for this to be worthwhile.

The error bands in Figure 12 represent the statistical variance of EVA set by the available observation periods, $O\approx$ 20 years, but which may contain some geographical dependence on climate. This is tested by mapping $f_{∆}$ in Figure 13, where: (a) for all stations shows some high and low values confined to, but interpolated around, individual stations, indicating probable statistical outliers; (b) culls all values outside the erf 5%−95% confidence limits, revealing a small east-west gradient.

4.1.5. Averaging Period Factor

The factor for averaging period, $f_{\overline{t}}$, can be directly estimated from the ratio of observations in strong winds over a given threshold, or else by comparing predicted extreme speeds at a given return period. Using the first approach, ensemble averages of all ASOS stations are listed in Table 5 for a range of threshold speeds and $\overline{t}=2^{\text{'}},10\text{'}$ and ${60}^{\text{'}}$. The population is very large for the lower thresholds, so there the statistical variance is small, and the standard deviations chiefly represent geographical variation. The population decreases rapidly above ~50kn, and for >60kn where the statistical variance is comparable to the geographic variance. Of most relevance is the factor to transform from hourly to 10-minute averages, $f_{\overline{t}}\{\overline{t}={60}^{\text{'}}\to {10}^{\text{'}},R=50\mathrm{y}\}$ and, using the second approach on the XIMIS estimates for $w=2$, the ensemble average was $〈f_{\overline{t}}〉=$ 1.117±0.0461 for all stations, reducing to $〈f_{\overline{t}}〉=$ 1.112±0.0368 when 5%/95% confidence outliers were removed, and this corresponds to a threshold speed between 50-60kn in Table 5. The [1] value $f_{\overline{t}}$ = 1.06, adopted into the Eurocode, has only half the increment above unity indicated by the ASOS data. Unpublished data for Marham, UK, indicates a lower value, $f_{\overline{t}}$ = 1.054, identical to [11]. This suggests that the synoptic scales of CONUS windstorms are about half that of the Atlantic depressions that dominate European extremes, and are compatible with the integral timescales in Figure 7 which are about half that found in the UK [28].

Figure 14 maps the station values of the factor, $f_{\overline{t}}\{R=50\mathrm{y}\}$ for $\overline{t}={60}^{\text{'}}\to {10}^{\text{'}}$: (a) for all stations and (b) for values within the erf 5%/95% confidence limits. This also reveals a small east-west gradient.

4.2. Rayleigh Simulation and EVA

4.2.1. Implementation

The Rayleigh 1000-year simulations enabled a detailed analysis of the parameters that governed the empirical ~20-year ASOS results using 50 times more data, leading to a substantial reduction in statistical variance. The fixed value of $w=2$ lies in the upper quartile of the observed ASOS range, and its suitability as the datum is confirmed later. Observation periods: $O=10$, 20, 50 and 100 years; integral timescales: $T=3$, 6, 12, 24 and 48 hours; sampling intervals $t_{∆}=1$, 10, 30, 60 and 180 minutes, and averaging periods, $\overline{t}=$ 2, 10, 15, 20, 30, 40, 60, 80 and 120 minutes, covered the ranges expected in practice.

4.2.2. Simulated Autocorrelations and Integral Timescales

The [35] and [17] methods produce instantaneous samples at 1-hour intervals when no averaging period is built into the respective target autocorrelation or spectrum. In a spectrum, the characteristic indicator of this is that the spectrum approaches the upper Nyquist frequency limit horizontally, since higher frequencies are aliased by reflection about this limit, as is evident in Figure 1 of [42]. In an autocorrelation, the characteristic indicator is its slope at zero lag – infinite (a cusp) when instantaneous, but zero (horizontal) when averaged. This is illustrated in Figure 15(a), showing the first $10\text{'}$ of lags for $\overline{t}=0\text{'},2^{\text{'}},{10}^{\text{'}}$ & $60\text{'}$ and target $T=12\mathrm{h}$. The decay of the autocorrelation moves to higher lags, raising the simulated $T$ by an increment somewhat smaller than $\overline{t}$, so the effect is more obvious at small $T$, Figure 15(b), than at large $T$, Figure 15(c). Averaging was achieved by simulating instantaneously at 1-minute intervals, then applying a running block-average filter with timescale $\overline{t}$. The simulated $T$ values for each combination of parameters are presented in Table 6. Note that the simulated $T$ values for $\overline{t}=0\text{'}$ match the target within 5% and the effect of $\overline{t}$ is of the same close order.

4.2.3. 50-Year Return Normalized Wind Speed

XIMIS was applied to each combination of $O,T,t_{∆},\overline{t}$ with $w=2$, producing multiple subsets of annual mode, $U$, and dispersion, $d$, as described in §3.5. For example, $O=10\mathrm{y}$ and $t_{∆}={180}^{\text{'}}$ produced 100 independent sets of 180 correlated subsets, yielding 18,000 estimates, whereas $O=100\mathrm{y}$ and $t_{∆}=1\text{'}$ produced only 10 independent single subsets. Ensemble averaging of the subsets ensured that the full 1,000-year record was employed in every combination to minimize analysis variance. As expected, dependence on observation period, $O$, proves to be minimal, as illustrated in Table 7 for contiguous 10-minute averages, $〈\overline{V}〉\{R=50\mathrm{y},\overline{t}=10’,t_{∆}=10’\}$. Here, it is seen that values for the four observation periods are identical up to the fourth significant figure over the range of integral length scales, $T$. This confirms that an ensemble of shorter observation periods is just as effective as one long observation period over the same record length, as previously suggested by [43] and endorsed by [44]. Consequently, values were further ensemble-averaged for all $O$, reducing the combination of parameters from four to three.

4.3. Rayleigh Disjoint Sampling Factor

4.3.1. Dependence on Disjoint Sampling Ratio and Timescale Ratio

Figure 16 presents examples of $f_{\mathsf{\Delta }}$ against $\mathsf{\Delta }$ evaluated for three values of timescale ratio, $\lambda$, and the datum return period, $R=50\mathrm{y}$, in the same format as Figure 12. Here, (a) and (c) show the single set of $T$ and $\overline{t}$ values that occur at the lowest and the highest evaluated $\lambda$. Both (a) and (c) are well outside the range encountered in practice but ensure that the practical range is covered. In the middle of the practical range, (b) for $\lambda =24$, contains four sets of $T$ and $\overline{t}$ that collapse together.

The fits to Equation (6), shown by the thick (red) curves, was made in two stages:

• Stage 1 – Optimization of coefficients $A$ and $B$ produced a stable trend with $B$ nearly constant for $6\le \lambda \le 72$, with 3 or 4 sets of $T$ and $\overline{t}$, but the single sets gave a high level of off-trend scatter.

• Stage 2 –Coefficient $A$ in Equation (6) controls scale whereas $B$ controls shape, so by fixing the value of $B$ at its median value in the stable range, $B=0.144$, the shape of the curve was fixed and $A$ was re-optimised.

• The model fit in (a) for $\lambda =1.5$ is poor, but this example represents a long averaging period with a very short integral timescale, implying ~2/3 of the windspeed fluctuations are filtered out, a combination which never occurs in practical observations. Conversely, all fits for $\lambda \ge 6$ are good, as in (b) and (c).

Figure 17 shows that the fitted values of $A\left\{\lambda \right\}$ collapse onto the power-law curve:

$$A\{\lambda ,R=50\mathrm{y}\}=0.3773{\lambda }^{-0.3184}$$

(7)

with a remarkable lack of scatter. This allowed the disjoint sampling factor, $f_{\mathsf{\Delta }}\{\mathsf{\Delta },\lambda ,R=50\mathrm{y}\}$, to be evaluated for any combination of $t_{\mathsf{\Delta }}$, $\overline{t}$ and $T$ through the $\mathsf{\Delta }$ and $\lambda$ ratios, as plotted over the practical field of observations on the contour chart, Figure 18.

Annotations A to H in Figure 18 represent observational parameters found in DS3505 and other national databases. For integral timescale $T=22\mathrm{h}$, typical of the UK and Europe: A = $10’$ means at ${60}^{\text{'}}$ intervals; B = $10’$ means at ${180}^{\text{'}}$ intervals; C = ${60}^{\text{'}}$ means at ${60}^{\text{'}}$ intervals; D = ${60}^{\text{'}}$ means at ${180}^{\text{'}}$ intervals. For $T=10\mathrm{h}$, typical of the USA ASOS observations: E = $10’$ means at ${60}^{\text{'}}$ intervals; F = $10’$ means at ${180}^{\text{'}}$ intervals; G = ${60}^{\text{'}}$ means at ${60}^{\text{'}}$ intervals; H = ${60}^{\text{'}}$ means at ${180}^{\text{'}}$ intervals. It is probably no coincidence that A and D lie on the same contour, each pair related exactly by ${\mathsf{\Delta }}_A/{\mathsf{\Delta }}_D=2$, ${\mathsf{\lambda }}_A/{\mathsf{\lambda }}_D=6$, and also for E and H, but the reason for this is presently obscure.

4.3.2. Dependence on Return Period

A benefit of adopting the median shape coefficient, $B$, in Equation (6) is then that $f_{\mathsf{\Delta }}$ for other values of $R$ are related to each other, and to the datum $R=50\mathrm{y}$, by a single scale factor. On plotting the excess over unity, $f_{\mathsf{\Delta }}-1$, in Figure 19(a) for all combinations of $\mathsf{\Delta }$, $\lambda$ and $R$ against the corresponding datum values, the results for each $R$ collapse onto a straight line through the origin, again with very little scatter. Values of the scale factor, defined by:

$$C=\left[f_{\mathsf{\Delta }}\left\{\mathsf{\Delta },\lambda ,R\right\}-1\right]/\left[f_{\mathsf{\Delta }}\left\{\mathsf{\Delta },\lambda ,R=50\mathrm{y}\right\}-1\right]$$

(8)

collapse in Figure 19(b) to an almost perfect onto a straight line when plotted against $\mathrm{l}\mathrm{n}\left(\ln \left(R\right)\right)$, to give the empirical fit:

$$C=1.355-0.2595\times \mathrm{l}\mathrm{n}\left(\ln \left(R\right)\right)$$

(9)

4.3.3. Analytical Model

Equations (6)-(9) combine to give a single empirical analytical expression for $f_{∆}$ which incorporates all the Rayleigh simulation parameters:

$$f_{∆}=1+0.3773\times {\lambda }^{-0.3184}\times \ln \left(1+0.144\times ∆\right)\times (1.355-0.2595\times \ln \left(\ln \left(R\right)\right))$$

(10)

Figure 20 tests the fidelity of Equation (10) against the simulation values from which it was derived. The small prediction error of ~0.4% is mostly from the statistical variance of the simulation, with a small contribution from the fit to Equation (7). There are >2,000 values in this plot, so the ~20 visible outliers represent only ~1%.

4.4. Rayleigh Averaging Period Factor

4.4.1. Derivation

Following the experience of the ASOS analysis in §4.1.5, the mean-over-threshold approach was discarded in favor of evaluation only from the Rayleigh simulation extremes, focusing on the factor $f_{\overline{t}}\{\overline{t}={60}^{\text{'}}\to {10}^{\text{'}}\}$ required to translate between the standard values: $\overline{t}={10}^{\text{'}},60\text{'}$, $t_{∆}={30}^{\text{'}},{60}^{\text{'}},180\text{'}$, at the datum 50y return period. As presented in Figure 21, transformation of the ordinate, $T$, into $\mathrm{l}\mathrm{n}\left(\mathrm{l}\mathrm{n}\right(T\left)\right)$, linearised the results as far as possible before fitting the quadratic expression:

$$f_{\overline{t}}=D+E\times \ln \left(\ln \left(T\right)\right)+F\times {\ln \left(\ln \left(T\right)\right)}^2$$

(11)

The fits are excellent, with unity correlations to the third decimal place. The coefficients D – F in Equation (11) are given in Table 8. This evaluation permits direct comparison with earlier studies. The two horizonal dashed lines represent the value $f_{\overline{t}}=1.06$ from [1] for Germany and $f_{\overline{t}}=1.054$from [11] for Denmark. The latter’s integral timescale values for 5 Danish sites, $〈T〉=19\mathrm{h}$ is comparable with [28] $T=22\mathrm{h}$ for Boscombe Down, UK, so that both these earlier estimates of $f_{\overline{t}}$ initially appear conservative against this study’s value, $f_{\overline{t}}=1.030$. However, this comparison lacks the additional factor needed to transform from hourly to continuous samples, $f_{∆}=1.049$, which brings the combined factors for $t_{∆}={60}^{\text{'}}$ to $f_{\overline{t}}\times f_{∆}=1.081$ for $T=22\mathrm{h}$.

The narrow focus of this evaluation is justified because $f_{\overline{t}}$ can alternatively be derived directly from the model for $f_{∆}$, Equation (10) over the full range of simulated parameters, due to the identities:

$$f_{\overline{t}}\{\overline{t}={60}^{\mathrm{\text{'}}}\to {10}^{\mathrm{\text{'}}}\}\equiv {\displaystyle \frac{\overline{V}\left\{\overline{t}=10\mathrm{\text{'}}\right\}}{\overline{V}\left\{\overline{t}=60\mathrm{\text{'}}\right\}}}\equiv {\displaystyle \frac{\overline{V}\left\{\overline{t}=10\mathrm{\text{'}}\right\}}{\overline{V}\left\{\overline{t}=0\right\}}}\times {\displaystyle \frac{\overline{V}\left\{\overline{t}=0\right\}}{\overline{V}\left\{\overline{t}=60\mathrm{\text{'}}\right\}}}\equiv {\displaystyle \frac{f_{∆}\left\{\overline{t}=10\mathrm{\text{'}}\right\}}{f_{∆}\left\{\overline{t}=60\mathrm{\text{'}}\right\}}}$$

(12)

for any given constant values of $t_{∆}$ and $R$. This is verified for the case of contiguous hourly-mean speeds, $t_{∆}=60\text{'}$, by Table 9 where the differences in value, which represent the overall model fitting variance, remain <1% except for the shortest timescale.

4.4.2. Contiguous Hourly Means to Continuous $10’$-Means

As $\overline{t}={10}^{\mathrm{\text{'}}}$ averages are the dominant datum for codes and regulations where convective gusts are not dominant, e.g. Europe [8], one principal need is to transform older, contiguous, ${60}^{\text{'}}$ averages into continuous ${10}^{\mathrm{\text{'}}}$ averages. It follows that transforming sampled hourly means to continuous $10’$ means requires two steps: a) for averaging period, $\overline{t}={60}^{\text{'}}\to {10}^{\text{'}}$; then b) for sampling interval $t_{∆}=60\text{'}\to 0^{\text{'}}$. The combined transformation factors for $R=50\mathrm{y}$ over the ranges of $t_{∆}$ and $T$ are given in Table 10, derived through Equation (12) by:

$$f_{\overline{t}}{f}_{∆}={f_{∆}}^2\left\{\lambda =T\mathrm{\text{'}}/10\mathrm{\text{'}},\mathsf{\Delta }=6\right\}/f_{∆}\left\{\lambda =T\mathrm{\text{'}}/60\mathrm{\text{'}},\mathsf{\Delta }=1\right\}$$

(13)

with $T$ evaluated in minutes, not hours.

From this it is seen that the previously proposed range, 1.054 – 1.06 for $T\approx 22\mathrm{h}$, underestimates by 2% – 2.5% through neglecting the effect of hourly sampling but are adequate for ½-hourly sampling.

5. Discussion

5.1. Variability of the Sampling Interval Factor

It is appropriate to comment on the lack of scatter in the fits shown in Figure 17 and Figure 19 and 21, which the casual reader might think is too good to be true. The asymptote of the upper tail of the FT1 extreme-value distribution is the exponential distribution, towards which exponential parents converge very, very quickly [29] (and conversely Gaussian parents converge very, very slowly). The XIMIS method exploits this rapid convergence by pre-conditioning the wind speed to be exponential, by transforming to $V^w$ in Equation (3), ensuring the most rapid possible rate of convergence to the asymptote. A defining characteristic of the asymptotic FT1 is that it maintains constant dispersion and shape, with only a linear shift to the mode in terms of $\mathrm{l}\mathrm{n}\left(R\right)$. The reduction in statistical variance in using a long, 1000-year, simulation is the same for any given $R$ whether this is treated as a single 1000-year record or as the ensemble average of ten 100-year records in this study. Taken together, these characteristics ensure very good fits with minimal scatter.

While the sampling interval factor$f_{∆}$, , scales by $∆=t_{∆}/\overline{t}$, as shown by [16], the ASOS observations and the Rayleigh simulation demonstrate that it also depends on the integral timescale – scaling by $\lambda =T/\overline{t}$, and also on the risk of exceedance – scaling by $R$. This means that a single value of $f_{∆}$ is only appropriate for regions where the wind climate and the integral scale, $T$, are effectively geographically constant. In nature, the wind climate is a mixture of meteorological processes, each with different timescales, so that the single-process Rayleigh model is appropriate at return-periods high enough to exclude the other processes. Where more than one process compete for the extreme at the relevant return period, typically $R=50\mathrm{y}$, the mixture must first be sorted into its individual components [45,46], when the model will apply to each component.

Equation (10) describes how the f-factors scale, but the values it provides are derived from the datum Rayleigh index, $w=2.$ This leaves the question of how sensitive these empirical coefficients are to other Weibull index values. As all the coefficients scale time while the index scales probability, the answer was anticipated to be “not by much”. This expectation is supported by the first-stage ASOS results which show identical trends and similar $f_{∆}$values within the error limits, but was confirmed by the additional sensitivity analysis reported in Appendix A. It would be instructive to repeat these calibrations using a simulation that accurately reproduces the seasonal-diurnal modulation of the mixture components of differing wind climates, e.g. as represented by the OEN mixture model [23].

5.2. Revisiting the ASOS Observations

Climatic conditions vary considerably across the CONUS as indicated by the major Köppen-Geiger climate classes [47], resulting in a mixture of synoptic-scale components acting in geographically-varying proportions, in addition to the convective gust and hurricane components excluded from this study. The major Köppen-Geiger climate classes shown in Figure 22 indicate that CONUS can be roughly divided into three major climatic zones:

1. West – Stations west of -100° longitude: principally Tundra and Desert classes, $〈T〉=9.0\mathrm{h}$, except along the Pacific coast;

2. Northeast – Stations east of -100° longitude and north of 38.5° latitude: Cold class, $〈T〉=10.5\mathrm{h}$;

3. Southeast – Stations east of -100° longitude and south of 38.5° latitude: Temperate class, $〈T〉=11.7\mathrm{h}$;

whereas $〈T〉=10.3\mathrm{h}$ for the mixture of all zones.

Figure 23 compares the ensemble-average $〈f_{∆}〉\left\{\overline{t}={10}^{\text{'}},{60}^{\text{'}}\right\}$ at ASOS stations for the datum DS3505 1-hour sampling interval against the Rayleigh model predictions in a Q-Q plot, over the range of $R$ for all stations and for each zone. In all cases of $\overline{t}={60}^{\text{'}}$, the values match well for $R=2\mathrm{y}$, but lie progressively below the 1:1 line as $R$ increases. This is consistent with the increasing influence of the EV dispersion parameter which is elevated by the higher statistical variance in the ASOS estimates at short observation periods, i.e., the cause of this discrepancy is probably in the observations rather than in the model. The same trend occurs in (a) for $\overline{t}={10}^{\text{'}}$ at all stations but shifted above the 1:1 line. Examining $\overline{t}={10}^{\text{'}}$ by zone, (b) shows a good match for $R=2\mathrm{y}$ but the values shift upwards progressively as the climate classes become warmer through (c) and (d). The principal difference between $\overline{t}={10}^{\text{'}}$ and $\overline{t}={60}^{\text{'}}$ is that the former admits additional mesoscale contributions, while the latter excludes them, and so their dominant mixture components are different. The progression from (b) to (d) is thought to reflect the increasing influence of thermally driven mesoscale components in the warmer climate zones that are not suppressed by the $\overline{t}={10}^{\text{'}}$ average. The integral timescales associated with these additional mesoscale components will be shorter than the values for the observed mixture.

By assuming the mesoscale components dominate the extremes for $\overline{t}={10}^{\text{'}}$, the corresponding values of $T\{\overline{t}={10}^{\text{'}}\}$ for the best fits to Equation (10) were estimated by optimization, as shown in Figure 24. The optimized fit for the western Tundra/Desert zone stations is only marginally different, indicating a good match with the observed $〈T〉$. The very good fits for the northeastern Cold and southeastern Temperate zone stations correspond to lower values of $T$, giving weight to the assumption that mesoscale components dominate in these zones. Applying this finding in practice would require prior knowledge of the value of $T$ corresponding to the dominant component when a mixture cannot be separated into its components, implying that a methodology to establish the dominant $T$ needs to be developed.

A conclusion can be drawn from this that the $\overline{t}={10}^{\text{'}}$ averaging period, as used in the wind code for Europe [8], is unsatisfactory in that it neither includes nor excludes mesoscale convective gusts completely. It seems sensible that 3s gusts should be used in locations where convective gust events dominate, as it is for CONUS [40]. Whereas hourly-means, with a gust factor approach, should be used where synoptic-scale storms dominate, as it was previously for the United Kingdom [7]. Like for CONUS, the climate types vary across the countries covered by the Eurocode [8] and, while this already allows for special zones in alpine regions, it would be sensible also to adopt a maximum gust approach to the thunderstorm-dominated regions of Italy [48].

6. Conclusions

• There can be no global constant value for the disjoint sampling interval factor, $f_{∆}$, owing to its dependence on climate, although a single value may be appropriate for regions where $T$ does not vary significantly.

• The variability of previous estimates is primarily due to climatic differences between study regions.

• The model for the disjoint sampling factor, $f_{∆}$, Equation (10), is appropriate when a single climate component is dominant, or when a mixture has been separated into its constituent components.

• Equation (10) is universal in the sense that, because it includes the effects of averaging time and probability of exceedance, the model for the averaging factor, $f_{\overline{t}}$, may be derived from it by Equation (12).

• The model for the combined averaging-sampling factor, $f_{∆}{f}_{\overline{t}}$, given by Equation (13), is for a Weibull index $=2$ (Rayleigh), which is expected to be appropriate for the dominant synoptic-scale component of a mixed wind climate.

• The sensitivity analysis (Appendix A) indicates that impact of index, $w$, over its observed range, is small.

• The model implies that the ten-minute mean and hourly-mean extremes in the ASOS observations for CONUS stations east and west of –100° longitude come from different mixture components, each characterized by a different integral timescale.

• A method to estimate the dominant timescale for extremes of a mixed wind climate needs to be developed.