1. Homogeneity of the Historical Time Series

13

Supplementary Material

1. Homogeneity of the historical time series.

First visual wave observations in the selected northern Black Sea ports were organized in 1893. However, regular routine visual observations in Yalta are available only since 1900, whereas in Odessa and “Chersonesus lighthouse” since 1915 and 1916 respectively. Until 1938 the 12-grade Beaufort scale had been used, and starting from 1938 it was replaced by 9-grade scale. The visual observations had been performed until 1954 when all maritime hydrometeorological stations were equipped by perspectometer-wavemeters. It should be emphasized that before 1913 the observations were performed by non-professional personnel. Only in 1913 the semi-annual courses for the observers were organized. So, in fact more or less qualified observations had begun in late 1913.

Taking into account the above-mentioned circumstances, the homogeneity of the time series must be carefully checked. In fact, there are at least two potential problems. First of all, initial non-professional observations can be principally biased and contain high-amplitude random noise. That is why the first 14 years of Yalta observations were excluded from analysis. Second, the Beaufort scale is nonlinear and discrete, and its use may introduce some bias in comparison with actual heights even if visual observations are performed by professional observers. Therefore we transformed the post-1954 actual wave heights into the Beaufort scale and tested the homogeneity of the time series for the pre- and post-1954 periods. We used 9-grade Beaufort scale for the entire period converting early 12-grade observations into 9-grade scale.

Figure S1 and Table S1 demonstrate the time series and some statistical properties of the stormy waves in Yalta in terms of the Beaufort scale. They clearly show that the routine observations’ data since the mid of the second decade of XX century characterize the wave heights and their variability quite reliably and they can be combined with the post-1954 semi-instrumental measurements. Small negative linear trend of the stormy wind-wave heights (Beaufort numbers) is not significant. At the same time superimposed multidecadal (~50 yrs) variability of the stormy wind-wave heights, seen in semi-instrumental post-1954 observations, also exists in the pre-1954 period. Thus, the low-frequency variability of the wave heights in the Northern Black Sea can be extracted using routine data since 1916, because they are quite reliable. In addition, data from all three stations are simultaneously available only since 1916.

It is remarkable that significant secular negative trend of monthly averaged wind speed and the associated decrease of monthly averaged wind wave heights were found on different Ukrainian stations (see e.g. (Ilyin 2009)). Thus our own and the published by other authors results show that the regional monthly averaged wind and Black Sea wind waves were weakening during the major part of XX century, in contrast to stormy wind waves. The heights of stormy waves are mostly characterized by quasi-periodical interdecadal variations because the linear trend of the maximum wind speed is much smaller than the linear trend of the monthly average wind speed (Ilyin 2009).

Table S1

Average Beaufort number and standard deviation of winter stormy waves for three periods of observations (showed also in Fig.S1) in Yalta in XX century

Yalta / Period / 1916-1946 / 1947-1971 / 1972-1995
Mean / 4.99 / 5.01 / 4.68
Sigma / 0.65 / 0.78 / 0.76

a

b

c

Figure S1. Variability of Beaufort numbers in winter of 1916-2010: Yalta (a), Chersonesus lighthouse (b) and Odessa (c). Lines show the decadal-scale trends for periods of averaging (see Table S1)

2. Statistics of extreme waves and selection of the stormy wave height threshold

Extreme wave height analysis is one of the useful tools to identify the stormy conditions. It uses the statistical functions of sampling waves including the highest waves (Hm exceeding definite threshold). However, the absolute scale of Hm seems to be not the perfect choice (Egozcue et al. 2005; Sanchez-Arcilla et al. 2008). Wave heights of ~7.0 m are similar to wave heights of ~7.5 m; in contrast, waves of 0.5 m and 1.0 m differ two-fold. In this case the log-transformation is a convenient scaling procedure. So subsequent analysis will be performed with a value X=log10 Hm as a measure of stormy wave magnitude.

There are few standard techniques enabling to carry out an extreme event analysis (Hawkes et al. 2008; Lopatoukhin et al. 2000; Van Gelder and Mai 2008). Two of them, viz. Annual Maximum Series (AMS) and peak-over-threshold (POT), will be used in extreme wave height analysis. The AMS method is based on the use of distribution of annual maxima which, for large samples, follows one of three asymptotic distributions: Gumbel, Weibull or Frechet. They are conveniently unified into a single Generalized Extreme Value (GEV) distribution:

,

where x is a shape parameter, s is a scale parameter and m denotes extremum localization.

In general case the AMS function depends on the type and parameters of the original wave height distribution. But if this type is unknown it is recommended to apply GEV (Coles 2001). A method of maximum likelihood (ML) is preferred to fit the distribution function. If the FAMS is defined, a return period of extreme log-wave-height can be written as

.

Parameters of GEV distribution estimated for the periods 1954-1983 and 1984-2010 are given in Table S2. According to the results of standard t-test there is a significant difference (on 95% significance level) between the parameters μ in the positive and negative WA phases for “Chersonesus lighthouse” and “Yalta”. Indeed, standard error for μ does not exceed 0.02. “Odessa” station is situated on the margin of wide shallow shelf. That is why the change in values of the statistical parameters between different phases is not so large.

Table S2

Parameters of GEV distribution function for different phases of WA

Station / 1954-1983 (positive WA) / 1984-2010 (negative WA)
m / s / ξ / m / s / ξ
Yalta / 0.55 / 0.10 / -0.31 / 0.40 / 0.08 / -0.87
Chersonesus lighthouse / 0.67 / 0.08 / -0.36 / 0.50 / 0.096 / -0.21
Odessa / 0.36 / 0.11 / -0.14 / 0.34 / 0.077 / -0.17

Note that the shape parameter ξ has negative values for both subperiods. This means that the GEV function can be reduced to one of its special cases – Weibull distribution function, which is characterized by a notable feature linked to the existence of an upper boundary of maximum wave height.

The AMS gives an understanding of extreme range of wave heights but is not suitable for analysis of storm events for long-term period. The reason is that AMS technique excludes the storms which are not the strongest ones in particular year but could be the strongest in other years (Lopatoukhin et al. 2000). This obstacle could be overcome by using POT method (Hawkes et al. 2008). Using the recommendation of the authors of the work Egozcue et al. (2005), storm event may be identified if observed wave heights exceed a certain critical threshold, u, for rather long period (say, 12 hrs or so for the Black Sea region).

In order to analyze the number of storm events (N) for different time intervals (say, t yrs) the statistical model of rare events (namely, the Poisson distribution) can be used (Egozcue et al. 2005; Sanchez-Arcilla et al. 2008):

, (1)

where λ is Poisson rate and is inversely related to the return period of events τ=1/λ. The stationary Poisson process is supposed to be formed by independent and identically distributed events. In order to test the difference between Poisson rates l1, l2 evaluated for the two periods one should use binomial test by conditioning on the total storm event count (Krishnamoorthy and Thomson 2004). In this case the ratio l1/l2 is a binomial success probability.

Each storm can be characterized by the maximum wave height (exceeding the definite threshold) being observed during certain period. Wave height is presented by log-transformed variable X. As a rule X is believed to have independent statistical properties from one event to another, but its distribution is supposed to be the same (Egozcue et al. 2005). Selection of the distribution function can be problematic but it is strongly recommended to use Generalized Pareto Distribution (GPD) applied to the variable X exceeded a reference threshold u or for Y=X-u (Hawkes et al. 2008). GPD function has the following form

,

where ξ is a shape parameter and β is a scale parameter (Egozcue et al. 2005). GPD function is flexible because it presents finite and infinite tailed distributions governed by two parameters. Three attraction domains (DA) of GPD function are recognized depending on the ξ parameter. If ξ=0 the GPD belongs to Gumbel DA and its support is infinite (ymax= +∞); if ξ>0 the GPD is in Frechet DA with infinite support (ymax= +∞); if ξ<0 the GPD is in Weibull DA and its support is finite with ymax=-β/ξ.

In order to specify the stormy conditions, the log-wave-heights must correspond to p values exceeding definite level. Usually the p value between 0.01 and 0.05 is being chosen. Two criteria of p selection are commonly applied. Firstly, the selected threshold must not be too high. Otherwise, the number of storms will be too low. The p-parameter corresponding to Kolmogorov-Smirnov goodness-of-fit test is widely used for these purposes (Lopatoukhin et al. 2000). Secondly, the log-wave-height limit must remove the low-energy events which cannot be considered as storms.

Other requirements are applied if GPD function is selected to estimate the extreme values. It is known that the mean excess over a threshold should be a linear function of storm intensity (Lopatoukhin et al. 2000) if shape parameter of GPD ξ<1,

.

So the reference threshold, u, should belong to “linear” segment of mean excess depending on the threshold values. In addition an inspection of “linear” segment allows us to check the fit of GPD and to define an attraction domain. Negative linear trend of the mean excess reflects Weibull DA, the positive one corresponds to Frechet DA. If the trend appears to be absent, hence the GPD is assumed to be in Gumbel DA.

The mentioned selection method can be applied to analyze the long-term variability of storm events in the coastal zone of the northern Black sea. In addition, one should keep in mind that storms can be observed on different stations simultaneously. Such cases will be considered as a single storm in the coastal zone of the northern Black Sea.

Inspection of Figure S2 using the above-listed recommendations clearly shows that the optimal range of thresholds begins from the values u>0.2 (Hm>1.5m). It is interesting to note that the value u=1.5 m is officially referred to as the threshold for dangerous event in Hydrometeorological Service of Ukraine. Slopes of the mean excess function are negative. This means that the Weibull DA is fit. The existence of the upper boundary of wave heights is physically valid in our case because of the semi-enclosed nature of the Black sea basin and dissipating processes acting in the coastal zone. In order to determine a threshold one should be guided by the goodness of fit of GPD. We put the reference thresholds for each station such as those listed in Table S3. These values define the number N of stormy wave events occurred during the period of 1954-2010.

Figure S2. Mean excess functions against threshold of log-wave-height to define a wave storm event. The right axis indicates the number of extreme events in the sample which is depicted as grey line.

Table S3

Parameters of GPD model for wave heights and the number of storm events (occurred during the period of 1954-2010) corresponding to them

Station / Yalta / Chersonesus
lighthouse / Odessa
u / 3.0 / 3.0 / 2.5
ξ / -0.305 / -0.34 / -0.16
β / 0.122 / 0.145 / 0.07
N / 115 / 341 / 62

Once the parameters of the GPD have been estimated, we could build a model of storm occurrence in time. As long as the storm event occurrence is assumed to be Poisson distributed (1), the probability of storm with a given intensity depends on the mean number, λ, of storms in a year. This is the key moment for the storm event evaluation. Poisson process is stationary which rules out account of any trends due to climatic variability in the datasets. However, wind wave conditions of the coastal zone of the northern Black sea are characterized by interdecadal variability. This means that the statistics of storm events may change drastically depending on the phases of WA process. This is the case for the return period which is defined as (Egozcue et al. 2005):

,

where τ(x)=1/ λ(x) is the return period of events with log-wave-heights X exceeding u.

3. Classification of large-scale atmospheric processes leading to the Black Sea storms.

One of the most effective methods to analyze macro-synoptic processes is their classification based on selection of common circulation patterns from a big variety of synoptic situations (Huth 2001). In this section, the classification of large-scale macro-synoptic processes for the storm events with wave heights ≥ 3 m (2.5 m for “Odessa” station) and 12 hrs life time threshold is discussed (for details, see section 2). It should be noted that the classification of large-scale atmospheric patterns for stormy conditions in the Black Sea has already been done in (Voskresenskaya et al. 2009), but without analysis of circulation types with respect to different WA phases. Here, we pay attention to large-scale (macro-synoptic) processes associated with storm activity in the Black Sea region and their occurrence in different WA phases.

The storm events with wave heights exceeding the chosen thresholds were observed most frequently on the “Chersonesus lighthouse” station. In total 341 storms were observed on this station in the period 1954-2008. 115 storm events were observed in that period on the “Yalta” station and 62 on “Odessa” station. Taking into account that storms may be observed on different stations simultaneously, the total number of storms in the coastal zone of the northern Black Sea, coming from all directions, in 1954-2008 was counted as 442. The procedure of hierarchical clusterization of 500 hPa geopotential topography patterns for these storm events permitted to extract four well-distinguished classes (or atmospheric circulation patterns). The circulation types were defined depending on localization of the mid-tropospheric ridge over the region restricted by 20ºN – 90ºN and 40ºW – 100ºE. The obtained principal spatial patterns are described below (for more details, see also (Belskaya 1949; Chernova 1959; Guijarro et al. 2006; Lepeshko 1989; Logvinov and Barabash 1982; Simonov and Altman 1991; Trigo et al. 2002; Voskresenskaya et al. 2009)).