Analysis of wave time series in the Estonian coastal sea in 2003–2014

Statistics of hindcast wave series were studied at eight differently exposed fetch-limited Estonian coastal sea locations (Harilaid, Kihnu, Kõiguste, Matsi, Neugrund, Sundgrund, Suurupi and Letipea). Based on episodic wave measurements obtained with the bottom mounted Recording Doppler Current Profiler in 2006–2014, a model for significant wave height was calibrated separately for those locations; using wind forcing data from Estonian coastal meteorological stations, a set of hourly hindcasts were obtained over the period from September 2003 to December 2014. Fourier analysis of monthly time series with regard to both average wave conditions and high wave events showed a distinct seasonality with specific amplitudes and the phases that differed by about a month between the locations. In addition, spectra of hourly time series indicated diurnal variations (24 and 12 h harmonics) and a consistent background noise following – 5/3 power law between the periods of 2 h and 4 days. Deseasonalized monthly data were used to analyse trends and extreme deviations (wavestorms). The possible trends, being obviously parts of longer variations, were not statistically significant over the studied 11 years. The lists of the most outstanding storms were different depending on the site.


INTRODUCTION
Statistical properties of coastal sea wave conditions are of great scientific and practical importance. Analysis of so-called wave climates usually includes a description of both typical (average) wave conditions and rough seas (Massel 2013). While tides are virtually absent in the relatively small and semi-enclosed Baltic Sea, waves, together with sea level variations during storms, become the major hydrodynamic agents (Soomere et al. 2008b) acting on seashores (Tõnisson et al. 2011) and shaping littoral habitats for marine life (Kovtun et al. 2011). The impact of windstorms and associated rough seas is frequently severe and can result in loss of life and extensive damage to property. Estonia lies in the zone of moderate latitudes (57.5°-59.5° N), which serves as a prolongation of the so-called North Atlantic storm track. Associated mainly with passages of cyclones, high wind and wave events occur almost regularly in every autumn or winter. In rare, extreme cases, significant wave height can reach about 10 m in the central Baltic Proper off the West-Estonian Archipelago (Soomere et al. 2008a). On the other hand, normal seasonality, as well as temporal variations in other scales, play an important role e.g. in hydrobiological processes of the sea, especially in boreal climates (Hünicke et al. 2015).
Wind waves have a certain amount of randomness. Subsequent waves differ in height, length and period with limited predictability even if their generation, propagation and transformation follow certain deter-ministic governing physics. Both wave models and measuring devices output their data as consecutive (e.g. hourly) values. In case of measurement, the various output values are already products of certain preprocessing and statistical analysis. The most widely used parameter is significant wave height (H s ). It is defined traditionally as the mean wave height of the highest third of the waves and nowadays it is determined through the wave spectrum. According to Rayleigh distribution, a one hundredth wave can still be 1.5 times larger than H s and one might even encounter a wave that is roughly twice the significant wave height (Massel 2013).
Climatological wave statistics can be computed either from direct measurements or wave model outputs. Both approaches have their strengths and limitations and have also undergone fast development during the last decade. Large numbers of wave modelling studies have been published on the northern section of the Baltic Sea (e.g. Soomere 2001;Jönsson et al. 2003;Alari & Raudsepp 2010;Räämet & Soomere 2010;Tuomi et al. 2011), and as a result, our knowledge of wave climates has certainly improved considerably. However, compared to a vast and growing pool of computer-generated hindcasts, only a tiny part has been properly validated against measurements. Although the commonly used models (e.g. Wave Predicting Model -WAM, Simulating Waves Nearshore -SWAN) have already proven themselves worldwide, a model together with a particular wind forcing at a particular application should be validated. That leads us to the importance of measurements, which basically can be either instrumental or visual.
The data obtained from instrumental measurements are seemingly 'correct', but even then a number of problems exist, e.g. due to cut-off frequency or possible malfunction of the equipment. Also, such data are essentially point-data, obtained from specific locations and do not easily provide spatial descriptions of wave conditions. The series tend to be short, much shorter than the time series available from weather stations or tide gauges. Continuous wave measurement is a demanding task and no long-term measurements (like the Swedish wave-buoy at Almagrundet, see, e.g., Broman et al. 2006) exist for the Estonian coastal sea. The data from currently operating four Finnish wave buoys (including one in the Gulf of Finland and another in the Northern Baltic Proper) are too sparse for Estonian applications (see Hünicke et al. 2015 for the most recent overview of the Baltic Sea wave measurements). As such devices should be removed before ice season (Pettersson et al. 2014), the measured wave data display extensive time gaps. Still, in recent years a number of bottom-mounted wave recorders have been installed near Estonian seaports (Alari et al. 2011), since typical and extreme wave regimes are of primary interest for navigation and harbour design. This means that they can mostly be found within harbours or marinas (e.g. Pärnu, Muuga, Rohuküla), where wave conditions are not representative for coastal or open seas. In addition, some episodic, experimental data sets exist.
We have performed a set of semi-empirical modelling studies based on wave measurements ( Fig. 1; Table 1) using a Recording Doppler Current Profiler (RDCP) oceanographic complex. The locally calibrated point model (LCPM) enabled us to produce multi-year wave hindcasts at a 1 h interval, as it would have been measured by the same RDCP at the respective locations. For the first time such a hindcast was presented for the Harilaid-Vilsandi region (1966-2006 by Suursaar & Kullas (2009). Hindcasts for Kunda-Letipea (Suursaar 2010) and some other regions followed (Suursaar 2013). While a shortcoming of the otherwise up-to-date spectral wave model application seems to be a somewhat limited reproduction of local wave properties in a shallow rugged coastal sea (e.g. Tuomi et al. 2012), largely due to relying on too general and 'smooth' gridded model winds, the LCPM applications may have their niche in wave studies (Suursaar et al. 2014), particularly in the shallow Estonian coastal zone where the swell component is usually small (Hünicke et al. 2015).
The first hindcasts using the LCPM were performed for the period from 1966 onwards, i.e., from the date when the digitized wind data were available, but it appears that such long hindcasts may be plagued by inhomogeneity of input data (Suursaar 2013;Suursaar et al. 2014). The exact influence of this inhomo-geneity is still unknown and therefore the present study analyses the supposedly highly homogeneous subset since September 2003. In addition to some previously discussed locations, Kihnu and Suurupi are entirely new sites included in this study (Fig. 2). The Sundgrund calibration is also new. All the existing hindcasts were updated until the end of 2014. In fact, they can be routinely updated as the calculation procedure is fast and straightforward for the pre-calibrated locations. We must mention, however, that no actual sea-ice cover is taken into account in this study. Consideration of ice in long-term wave statistics poses some difficulties (e.g. Tuomi et al. 2011). One should decide whether to assign zero values or omit such days, and then the question of how to treat the correspondingly shortened time periods in statistical analyses arises. Hence, the wave hindcasts here represent 'climatological' wave conditions, regardless of ice. Indeed, while some of the locations (e.g. Kihnu, Matsi and Letipea) may have been covered by ice for up to several months, mostly from January to March, others see little or no ice during most of the winters. Moreover, due to a significant increase in winter air temperatures over the last fifty years, the ice conditions in the Baltic Sea have gradually turned milder as well (Leppäranta & Myrberg 2009). Typical duration of ice cover has shortened by about two months in the Gulf of Riga and by about a month in the Gulf of Finland (Jaagus 2006), and the trend is expected to continue (Luomaranta et al. 2014). Among the studied eleven years, six winters were mild, including the two extremely mild winters of 2013/2014 and 2014/2015, when notable ice was not encountered in either of our study areas.
The objectives of the study are to (1) discuss our experiences with locally calibrated wave hindcasts in eight differently exposed coastal sea locations; (2) update the previously published hindcasts until the end of 2014 and to introduce three new ones; (3) perform time series analysis of hourly hindcasts from September 2003 to December 2014, including isolation of regular cycles, identification of peaks and examining possible trends and (4) discuss the climatological background of the obtained findings in relation to different expositions and local features.

Semi-empirical hindcast of waves
The study stems from long-term wave hindcasts performed for the selected coastal sites using a semiempirical model. The model was independently calibrated against the earlier wave measurements in those locations, thus, the first important pre-requisite was the availability of suitable wave measurements. Since 2003, hydrodynamic measurements using the RDCP-600 oceanographic measuring complex (manufactured by Aanderaa) have been carried out in 12 locations over the Estonian costal sea during 20 measuring campaigns with a total duration of about 1600 days. Eight locations included the measurements that have been suitably placed for nearshore wave conditions and also were extensive enough (i.e. one cut at least about 40 days long) for wave model calibration. The sites were all 1-3 km off the nearest coast and the mooring depth varied between 10 and 20 m. This is the depth range where the bottom mounted instrument, including its high-accuracy pressure sensor, can yield reasonably good wave data. On shallows, the waves are restricted by small depth and distorted by coastal effects. At deeper locations, a considerable part of the wave spectrum becomes damped (cut off) and altogether different instruments (e.g. surface buoys, wave-riders) should be used. Although the damping of the dynamic pressure with depth is compensated by the RDCP software, the quality of such measurements may deteriorate if the depth becomes larger than about 20 m. Wave-riders (normally used in the open sea) and pressure-based instruments (used in coastal and harbour applications) apply different principles for measuring waves and it is not possible to decide which way is more accurate. Their scope of use is somewhat different. The RDCP output preciseness is nominally 1 cm for wave heights and 0.01 s for wave periods. The 1 h measuring interval was set the same as it was in the  Table 1). The locations of the weather stations used in calibrations or hindcasts are marked with . routinely measured meteorological data used in the calibration and hindcast. More details on the measurements are available in Suursaar et al. (2012Suursaar et al. ( , 2013. A simple, fetch-based wave model was chosen for hindcast. Also known as the significant wave method, presented e.g. in the Shore Protection Manuals of the U.S. Army Corps of Engineers (USACE 1984), it calculates the significant wave height, period and wavelength as a function of wind speed, effective fetch length and water depth. The specific formulas used in this study were exactly the same as presented in Kullas (2009) andSuursaar (2013).
The fetch is usually measured as the headwind distance from the nearest shore for the applicable wind direction. In fact, the manuals and handbooks include a wide choice of such equations with slightly different empirical coefficients and procedures, which are supposed to take different basin properties and wind conditions into account (Bishop et al. 1992;Massel 2013). Some authors have developed GIS-based procedures for determining openness (fetches) for difficult coastlines (e.g. Tolvanen & Suominen 2005). In practice, it is still difficult to match the exact influence of islands, shoals and the coastline on actual waves. Our intention was to calibrate a wave model using high-quality wave measurements, so that afterwards it can reproduce the wave series for the period when the instrument was not moored to the sea but wind forcing data were still available from the source station where the model calibration originated. We admit that the choice of the exact formula version from the family of fetch-based models was not crucially important. Instead, the sitedependent calibration procedure became decisive for the model set-up.
For supplying the wave model with wind speed and direction, we acquired data from the meteorological Correlation coefficients (r) and rootmean-square errors (RMSE) between the corresponding raw data series are given in legends. Note that the measured/modelled peak on 13 December extends for 5.1/4.6 m (C). stations operated by the Estonian Environment Agency (EEA). For each of the wave modelling sites, we used wind data from the closest (coastal) station as the first choice. Data from four stations were used in this study (Table 1). While Vilsandi and Kihnu stations successfully reproduce marine or semi-marine wind conditions, the stations at Kunda and Pakri offer somewhat mixed quality in that sense (Keevallik et al. 2007), but are still the best available wind sources for the nearby wave applications.
The model calibration included a search for the appropriate depth parameter. This kind of a simple model assumes that the basin has a constant, generalized depth. The water depth in the model should thus represent both the depth of the actual mooring (i.e. 10-20 m) and the average depth of the nearby subbasin. Our calibrations yielded depths between 19 m (Kõiguste) and 33 m (Neugrund) (Suursaar 2013). Still, the most influential item in the calibration procedure was prescribing the fetch for different wind directions. After initial (rough) measurement of fetches from nautical charts, a series of new distributions of fetches was eventually created by maximizing the correlation coefficient and minimizing the root-mean-square error (RMSE) by means of consecutively adjusting the fetch in all 20° wide sectors. By trying to keep the maximum and average wave heights equal in the modelled and reference series, the iterative calibration procedure found the best set of fetches and, in a way, compensated for local wind impediments around the specific weather station. Wind forcing data derived from the coastal measurements are usually far from ideal, which should mean 'open terrain' or full openness to every direction. Calibration provided feedback on how the specific wind forcing might yield the actually measured waves and modified slightly the geographical fetches correspondingly. As just one example, Pakri wind data projected at the Neugrund location required the enhancement of wind speed (i.e. fetch) for the SSW direction. The direction was somewhat shaded by local vegetation and buildings around the Pakri meteorological station and eventually required up to 200 km long manipulated fetches to act properly on Neugrund waves.
Calibration quality varied depending on the complexity of the bathymetry and coastline, openness of the terrain around the wind instrument, but also on the distance between the wave and wind measuring locations. The prognostic value of calibration was higher when the calibration period included a wide range of wind directions and speeds. In fact, in most cases the measurement was extensive enough and even included some severe winter storms. Still, calibration efforts at some locations (e.g. Sillamäe, Küdema) were never precisely developed into long-term hindcasts because the measurement period did not provide such a variety in wind and wave conditions. In best cases, the correlation coefficient between the unsmoothed hourly data exceeded 0.9 (at Harilaid, Kihnu, Letipea, Kõiguste and Matsi) and the RMSE was around 0.2 m. We stress again that the mean and maximum heights of measured and modelled series were kept equal and the comparison was performed between statistically non-manipulated raw data (see Fig. 2).
The calibration details for some of the locations have previously been presented by Suursaar (2010Suursaar ( , 2013 and Suursaar et al. (2012Suursaar et al. ( , 2013. The new, Kihnu calibration, was quite successful (Table 1; Fig. 2A), yielding a very good correlation and small bias. The Suurupi calibration was for some reason much more difficult. It yielded moderately good comparison statistics (Fig. 2C), but most importantly, it was not possible to correctly reproduce the measured 5.1 m high H s peak during a storm on 13 December 2013. The model gave just about 4.6 m, indicating that certain nearshore phenomena (breaking, interference or backscatter) could have occurred during this prominent storm, requiring perhaps a special study. The Sundgrund calibration (Fig. 2B), where Pakri data were used instead of winds from the Lääne-Nigula station , is also new. Although Pakri is probably the best available long-term coastal station in NW Estonia, its winds may be plagued by the influence of the North Estonian Klint (Keevallik & Soomere 2009). The location of the station on the Pakri Peninsula has changed three times, but this fact does not interfere with the current study which covers the period since its last change in September 2003.
In addition to generally good calibration results, a few independent validation examples between the precalibrated model runs and field measurements can be found (e.g. fig. 3b in Suursaar 2010). An additional comparison (validation) between the outputs of the LCPM and the SWAN model forced with BaltAn65+ reanalysis (Luhamaa et al. 2011) wind data was presented by Suursaar et al. (2014). The overlapping time period of two independent wave hindcasts  enabled the models to be cross-validated in three different timescales: hours-days, months and years. The LCPM output retained a very similar visual appearance and statistical structure to the RDCP measurements, both having a 1 h time step. Owing to the 3 h time step and inherited temporal properties from the BaltAn65+ reanalysis wind data, the SWAN output was smoother and less detailed, but the difference (bias) between the models was small. The models showed a good agreement both in shortterm and monthly scales, with the exception of the multi-annual scale, where inhomogeneity issues of input data start to play a considerable role. It seems that the long-term wind input for the point-model had a small decreasing, artificial trend component. Although we tried to eliminate the influence of the most obvious instrument change of 1976 in the older data, there may still be other smaller inhomogeneity sources left. Starting from September 2003, the input data should be homogeneous over time, except in Kunda. As the wind instrument was transferred to the outer pier of the Kunda port in March 2014, which is windier than the previous location, the homogeneity of that series was compromised after that.

The study locations and the time series
The eight locations have somewhat different bathymetries, expositions and fetch distributions, which means that wave growth is differently selective to wind conditions in those areas ( Fig. 1; Table 1). In medium-sized water bodies, larger waves are favoured to propagate from the direction of the longest fetches. From the direction of a very short fetch, no matter how strong the wind can blow, wave growth beyond a certain (very low) limit is not possible (Massel 2013). Bottom topography in most of the 10-20 m deep study locations, 1-3 km off the nearest shore, was gently sloping towards the longest fetches. Only at Neugrund, the 100 m isobath lies just about 7 km away. This fact, together with the rugged bottom topography of the Neugrund submarine impact structure (see e.g. Suursaar et al. 2013), somewhat reduced the quality of wave modelling.
One should bear in mind that based on specific hindcasts, one cannot automatically decide which region generally has higher waves. Strictly speaking, the hindcasts were representative for the exact locations where the wave measurements and calibrations were initially performed. Moving away by just a few kilometres (e.g. to a deeper location) may significantly alter the wave height statistics. However, we assume that the hindcast results are suitable for studying long-term changes in wave conditions, which in turn should be representative for a wider area with similar exposure.
Technically, all the time series could start from 1 January 1966, as from this date, the digitized wind data are available in the database of the EEA. However, it is known that the wind measurement procedure and instruments have changed several times since then (see Keevallik et al. 2007 for an overview). Shortly put, winds were measured by 'weathercocks' (wind vanes of Wild's design) during 1966-1976 and automatic anemorhumbographs during 1976-2003, while MILOS-520 automatic weather stations have been operating since September 2003. The data from January 1966 to August 2003 have a time interval of 3 h and they are also less exact. The currently used MILOS-520 provide a 0.1 m s -1 value step and 1° directional resolution output in two versions. The 1 h averaged hourly wind data were used in modelling, as they yielded marginally better wave calibration results than the 10 min average data recorded once an hour. In general, the difference between these two wind instrument outputs is very small. We did not use wind gust data in modelling.
Old instruments were replaced by automatic weather stations in Estonia in September 2003, which still yielded about 99 000 hourly recordings until the end of 2014. The completeness of the data series is very good, reaching 99.9% in most cases. We used the hourly data in spectral analysis and in exploratory analyses, including distribution fitting (Fig. 3). In some cases, it was still more convenient to operate with monthly time series, which can be obtained starting from October 2003 (or January 2004, if complete years were preferred). For statistical purposes, we calculated summary statistics (e.g. average, standard deviation, 95th percentile, 99th percentile and maximum H s ) both on an annual and monthly basis. As it is common for wave models, the statistics were all based on significant wave height measurements (H s ), which are still conveniently close to the visually perceived 'average wave heights' (Soomere et al. 2008b). We also admit that it might be difficult to comprehend the true spectral meaning of e.g. the 99th percentile of H s samples. When creating monthly time series from actual hourly measurements for Fourier analysis, the length of a month was taken as 30.3 days with the central point of each element at the beginning of the 16th date. For descriptive statistics, trend analysis and Fourier analysis, the STATISTICA software package was used.

Seasonality
Time series formed on the basis of monthly statistics (averages, 95th and 99th percentiles) included n = 132 elements and covered 11 full years from January 2004 to December 2014 (Fig. 4). The only exception was Letipea, where due to the change in wind measurement location, 10 full years until December 2013 were considered. Depending on the calendar month, each element represented averages of 672 to 744 hindcast values. Thus, calculated on the basis of H s , the 99th percentile represents 7 h with the highest waves in a month. Single maxima could still be much higher, but carrying a considerable random imprint, such maxima were not included in the analysis of seasonality. In general, monthly 95th percentiles exceeded monthly averages by 2.6 times and 99th percentiles by about 3.5 times (Table 2).
Besides the possible inter-annual variability and trend, the most distinct feature of all series was seasonality. Peaks appeared more or less regularly in every autumn or winter and lows in early summers. Using Fourier analysis of monthly series, we can pick up the component describing seasonality and compare its amplitude, importance (i.e. share in total variance) and phase. We assumed that the seasonal cycle is not quite similar in different locations.
Time series can be represented as a sum of 2 n simple sine waves (harmonics): where 0 2 A is an average of the time series, i A and i B are the amplitudes of harmonics ( 1, , 2), i n   t is time (in months) and it represents argument function (angle in radians). STATISTICA software provided A and , B but for more convenient form 2 0 the amplitude i C and phase shift i  can be computed as follows: 2 2 ; arctan .
The first Fourier component had the same period length as the series (132 months), the 66th had a period of 2 months and the 11th was the one with the period of 24 h (10th in case of Letipea).
Seasonality in average wave conditions seems to be most pronounced at the locations surrounding the semienclosed and medium-sized Gulf of Riga. It is reflected both by handsome amplitudes and substantial shares of such motions in total variances. The amplitude of seasonal harmonic (C, Table 2; Fig. 5) was typically about 0.15 m (or up to 0.3 m for the annual ranges). Due to much higher values, the amplitudes were also much higher for the high percentiles. The seasonal harmonic explained about 40-70% of the total variance in Kõiguste, Kihnu and Matsi (being maximal 74% for the 95th percentile at Kõiguste). The share was somewhat smaller (10-40%) at the locations facing North, West and the Baltic Proper (Table 2). However, even then the seasonality remained a prominent and statistically significant feature of the time series, as its share well exceeded 1.7% (i.e. 1/60) expected from the white noise spectrum. There can be some discrepancies between the empirical seasonality and its idealized, sinusoidal form (e.g. Fig. 5B). Despite some idealization of seasonal processes, the few constants found from harmonic analysis can be useful for clarifying certain important, but otherwise not so obvious properties of the time series, in reconstructing time series and exclusion of regular cycles in trend analysis.
Seasonality did not follow the same timings; the phase in average wave conditions was shifted by up to 1.5 months in different locations and up to about a month in high wave conditions. On the basis of ,  we calculated the exact statistical dates for maxima and minima, assuming that seasonality was a sinusoidal process ( Fig. 5; Table 2). On a monthly basis, the average wave heights were usually the highest in December and the most common month of maxima for 95th and 99th percentiles was November. The empirical minima can Fig. 3. Empirical distribution function of hindcast wave heights at Harilaid together with a theoretical distribution curve. While short-term distribution of different single wave heights can be described using Rayleigh distribution, the long-term bulk data usually satisfactorily follow either log-normal or Weibull distributions.  Table 1 for codes of the stations. occur irregularly between April and July with specific dates in May or June, according to the corresponding sine components ( Table 2). The peaks were early at Matsi and Kihnu and the latest at Letipea. The reasons behind the specific expressions of seasonality in wave conditions including the different phases were climatological. Although selectively with regard to different wind directions, the local wave conditions always somehow reflect the corresponding wind climate (Räämet & Soomere 2010;Suursaar 2013).
Estonian climate and its temporal variations are largely controlled by the balance between the two competing (i.e. North Atlantic Oscillation -NAOdriven zonal and meridional) climatological phenomena (Jaagus & Suursaar 2013;Suursaar et al. 2015). Higher cyclonic activity and a more maritime climate dominate at Harilaid, Sundgrund, but also at Kihnu. The northerly exposed and easternmost Letipea is under the influence of anticyclonic, more continental climates. Indeed, the seasonal wind roses are spatially rather different in Estonia (Jaagus & Kull 2011). Moreover, directional distributions of strong winds substantially differ from mean wind conditions. At coastal stations, strong winds tend to blow from the directions of vast spans of the open sea, as the friction above the sea surface is considerably lower than above the land surface. Such statistics are not constant in time. According to Jaagus & Kull (2011), the percentages of westerly and southwesterly winds have had clear positive trends in Estonia over the last decades -especially in winter. For instance, at Vilsandi the directional winter-time peak was southeast in the 1960s, but southwest in the 2000s.  A, B, C) found from Fourier analysis for monthly averages (Av), 95th and 99th percentiles. 2C would be the range for seasonal variations. The calculation of dates for annual maxima and minima was based on the phase (), considering 15 January as the initial point of the series and the average length of month 30.3 days. S -share of the seasonal component out of 66 Fourier coefficients (with total variance 100%). * Calculation based on 10 years and 60 harmonics. Location codes explained in Table 1 Location code

Spectra of hourly time series
Fourier analysis of hourly time series (n  99 000, depending on the exact beginning date of the series in September 2003) yielded in total 2 n coefficients, provided that no time-series padding was used. The padding option is sometimes used when data are not collected at exactly regular intervals or include gaps. Such series are called 'dirty' and padding essentially increases the number of coefficients. In our case, the series were nearly perfect. For instance, 49 620 harmonics were obtained at Kõiguste (Fig. 6A). The sampling rate 1 cph (cycles per hour) gives the frequency for the shortest harmonic of 0.5 cph or 0.0018 Hz (corresponding period 2 h). The spectrum can be visualized both as a function of frequency or period. The simple spectrum (periodogram) obtained from Fourier analysis is a very noisy function (Fig. 6A). By averaging over a number of periodogram values (using Hamming weights in our case), a smoother energy spectrum can be obtained ( Fig. 6B-D). Because the number of coefficients was huge and both the frequency and variance span over several orders of magnitude, logarithmic presentation of the axes was convenient. We also normalized the spectra in such a way that total energy (variance) was equal to 100% in each case. Although the total variance differed up to threefold at our study locations, the normalization made the spectra directly comparable to each other.
The spectra showed some peaks over background noise which followed the so-called power -5/3 law (i.e., the function appears as a line with that slope when logarithmic scales are applied) in the band with periods between 2 h and about 4 days (Fig. 6A). Depending on the scale of the motions and the specific processes involved, the slope of the power law can be different. Between 4 days and one year, it was -1/3 and over the whole spectrum it was likely -1 (pink noise). Indeed, a family of pink-like and red noises (which also includes the -5/3 noise) widely occurs in nature (e.g. Ozmidov 1965); such noise is even considered as nearly ubiquitous. Red (Brownian) noise is considered as appropriate background noise for many climatic signals. The -5/3 background spectrum, for instance, also appeared in the band between 20 min and 4 days of the flow measurements in the Suur Väin Strait (Suursaar et al. 1995). The somewhat elevated high-frequency band (periods between 2 and 4 h; Fig. 6A) included the effect of aliasing. It is a common artefact in such time series where the motions shorter than the Nyquist frequency create false frequency, which is added up to the nearby 'visible' spectrum. Carrying very little energy, the highfrequency band is not particularly interesting in our case. Actually, half of the harmonics, a huge number, shown in green in Fig. 6A, comprised just 1.3% of the total energy.
Besides the seasonal peak analysed in the previous chapter, a diurnal cycle with the main period of 24 h was clearly visible in all the time series. Relatively, its proportion was highest at Suurupi, followed by Kõiguste and Harilaid. It originates from the diurnal variations in wind conditions (e.g. Keevallik & Soomere 2009), and most notably from the switching of breeze directions. Occurring mainly during the warm half-year, the breeze is a local phenomenon and it manifests itself better where larger thermal contrasts between land and sea can build. Although breeze-generated waves can occasionally be found to a distance of up to 100 km from the coast (Vethamony et al. 2011), it is mainly limited to quite a narrow coastal sea region (e.g. Soomere et al. 2012). It is also distinctive that the 12 h harmonic originated from the main diurnal cycle, while the further respective harmonics (6, 4 and 3 h) did not differentiate from the background spectra.
The important peak representing seasonal variations had the largest share at Kõiguste, where the corresponding single (T = 8760 h) harmonic of non-smoothed spectra (Fig. 6A) comprised as much as 9% of the total energy. Seasonality was relatively weak (but still visible) at the northerly exposed Suurupi location (Fig. 6D). On the other hand, the Suurupi spectrum displayed quite high energy in the low-frequency band. We can also notice a somewhat larger share of motions within the band of 20 days to 6 months at Harilaid, and a relatively lower share of such motions in the other locations.
Finally, by examining the spectra we are convinced that our hindcast data truly reproduce certain fundamental geophysical processes. The data were produced by a semi-empirical model which, indeed, faithfully mimics the behaviour of actual wave measurements. It is also a proof for the hindcast quality -as long as the wind input data itself are homogeneous over time. It seems to be a relatively recent understanding that the largest source of uncertainty in contemporary wave modelling is actually the wind information (Hünicke et al. 2015). A correct model alone cannot yield correct results if the input data are erroneous or nonrepresentative.

Trends
Depending on initial points and data treatments, linear trend analysis showed slightly different trend values for different time series versions (Table 3). For average wave heights, the annual changes by trend varied between + 0.12 and -0.73 cm, yielding + 1.2 to -7.3 cm total changes for 11 years. The changes for high wave conditions were larger (up to about 1.5 cm per year in either direction), but still not quite high, as the full variability range reaches up to 5 m (Fig. 4).
It should be noted that trend analysis of monthly data should be generally avoided, even if full years are covered (e.g. 132 months from January 2004 to December 2014). For instance, if according to typical seasonal cycle the highest values should be expected in December or November (Fig. 5), the series which begins with an element from January and ends with December may include an artificial trend component due to difference between the expectances of December and January values. However, knowing Fourier coefficients for seasonality, it is possible to subtract this regular component and analyse the residuals. The procedure lowered the variance by about a half, but it did not eliminate higher harmonics of seasonality or meteorological variability. It is also possible to analyse the series which have not exactly full year's lengths.
The lengthened 135 month versions that started from October 2003 did not introduce any significant changes to trend courses. In fact, none of the trend slopes (Table 3) were statistically different from zero (on p < 0.05 level). Although some trends seemingly did exist, the monthly and interannual variability of series was simply too large (Fig. 4). Also, the 11 year period is climatologically too short for showing clear tendencies and it can be viewed as a part of longer variations (e.g. 20-30 year quasi-cycles; Jaagus & Suursaar 2013). According to many climatological studies, regime shifts associated with changes in cyclonic activity in the region of the Baltic Sea (Stips & Lilover 2012) occurred in 1989 and the mid-2000s (down). Corresponding high phases in the 1990s-2000s can be followed both in sea level and wave statistics , as well as in the intensity of coastal change (Tõnisson et al. 2011).  Table 2 for the exact parameters and Table 1 for codes of the stations. Fig. 6. The normalized energy spectra at Kõiguste (A, periodogram, a function from frequency ω; B, energy density spectrum, a function from period T); energy density spectra at Harilaid (C) and Suurupi (D) together with some characteristic peaks and slopes of background spectra. The red circle (in A) indicates the influence of aliasing. Table 3. Trends (cm year -1 ) in monthly data of averages (different versions), 95th and 99th percentiles. Monthly data from January 2004-December 2014 or October 2003-December 2014 (*deseasonalized); annual averages include 11 years (** 10 years and correspondingly 120 or 123 months in case of Letipea). Location codes explained in Table 1 Annual averages Monthly averages Monthly 95% Monthly 99% Location code 2004-2014 2004-2014 2004-2014* 2003-2014* 2004-2014* 2004- The prevalence of declining trends over the last 11 years (Table 3) may indicate that the NAO-driven cyclonic activity and the westerlies are currently in the phase of reduction above Estonia. On the other hand, the trends for high wave events seem to be increasing at northerly exposed Letipea (Fig. 7). When westerlies increased in the 1990s and the corresponding wave heights rose at Harilaid, the northerly winds decreased at Kunda and the waves considerably decreased at Letipea (Suursaar 2010). The fact that changes at the westerly and northerly exposed locations are in contraphase could be explained by the influence of long-term variations (i.e. north-or southward shifts) in typical cyclone tracks above the Baltic Sea (e.g. Pinto et al. 2007). There were more cyclones in the 1990s, which bypassed Estonia from the north, creating strong westerly winds and less northerly winds. Besides, in decades when there were more such cyclones that crossed Estonia, increase in northerly and decrease in westerly winds could be observed. In terms of annual resulting wind directions (calculated from the wind velocity components), the prevailing wind directions have changed from about 220° to 230° at Kihnu, from 230° to 240° at Vilsandi and from 190° to 210° at Kunda in 1966(Suursaar 2013.

High wave events
Wavestorms associated with high wind events were analysed on the basis of monthly 99th percentiles as well as in individual maxima from the wave hindcasts. The 99th percentile in statistical terms reliably describes a sufficiently long (at least 7 h) and influential storm. Although a single maximum value can definitely be found within a violent storm, its exact numeric value is not very reliable for mainly two reasons. Firstly, model performance can be somewhat inexact in case of very high winds (and waves) extending beyond the variability range which was available during the model calibration. Secondly, especially during violent storms, occasional disruptions of the wind measuring equipment are quite frequent and the actual wind speed maxima may have not been recorded. Wavestorms typically occur in autumn or winter, which are windy seasons in the Baltic Sea region. Actually, some kind of high wave conditions could be statistically expected on an annual basis ( Fig. 4; Broman et al. 2006). In search for extreme deviations from the 'norm' we subtracted the corresponding regular seasonal cycles from the 99th percentile series and identified the most outstanding wavestorms for each location. Even then, larger deviations mostly occurred during the cold half-year. Not only are waves higher, but also variability is higher in autumn and winter months. The randomness of high wave events was higher in northerly exposed coasts (Letipea) where the longest fetches coincide with the northern directions where strong winds are infrequent. On the other hand, westerly stormwinds on the west coast can be expected every winter.
Although a few major storms can be found in the records for most locations (e.g. a storm on 9 January 2005), the lists were still sufficiently different (Table 4). A curious exception was the southeasterly exposed Kõiguste location, where none of the regionally important storms left a mark. Instead, it was sensitive to rare and not so influential easterly wind events which were even less influential for westerly exposed locations. Temporal variations of positive residuals (Fig. 7) showed several subtypes, evidently owing to different local exposures. At westerly exposed Harilaid, but also at Kihnu, a decrease in high wave heights has occurred. At northwesterly exposed Sundgrund (but also at Neugrund and Suurupi), no clear linear trend was observed. Instead, certain subcycles with high phases can be found around 2004 and 2009 and a low phase in between. At the easternmost and northerly exposed Letipea, the violent northerly wind-storms seem to occur irregularly. Their strength may be quite unexpected and the corresponding year-to-year variability of wave conditions is huge. In a way, at the northern coast, regionally dominating and NAO-driven westerly activity mingles with the locally favoured effects of northerly exposure . In addition to the locations represented in Fig. 7, the months with some highest waves were as follows: The single maxima mostly occurred during the above-mentioned months as well, but not necessarily in the same order (Table 4). In addition, a few more wavestorms should be highlighted. A modelled 4.67 m value (from 13 December 2013) from Suurupi shortly falls off Table 4, but the actually measured H s value of 5.11 m (see also Fig. 2) seems to be one of the highest significant wave heights actually ever measured in the Gulf of Finland (e.g. Pettersson et al. 2014). In addition to high H s , a 7.3 m maximum wave height and a peak wave period of 10.9 s were recorded using the RDCP. It is still hardly reproducible using the given winds (storm Billie; 21.8 m s -1 hourly average westerly wind at Pakri, gusts up to 31 m s -1 ) and the 20 m deep location. Another value not included in the analysis is a 19 m s -1 Table 4. Values and dates for two single highest (H s ) modelled wave events for each location together with the corresponding hourly sustained (and gust) wind speed (m s -1 ) and peak wind direction. Location codes explained in Table 1 Location code

CONCLUSIONS
Statistical analysis indicated that the semi-empirical model together with its fetch-based calibration scheme can deliver reasonable wave hindcast results which reliably mimic the behaviour of actual wave measurements in those locations. A set of high-quality hourly hindcasts obtained at eight differently exposed locations over the period from September 2003 to December 2014 enabled the study of statistical properties (trends, spectra, peaks and regular cycles, such as seasonality and diurnal variability) for both mean and high wave conditions. All the existing eight long-term hindcasts can be routinely updated in the future, provided that the wind input remains homogeneous in a statistical sense. Fourier analysis of time series comprising monthly samples with regard to both average wave conditions and high wave events showed a distinct seasonality with specific amplitudes and the phases that differed by about a month between the locations. In addition, spectra of hourly time series indicated diurnal variations (24 and 12 h harmonics) and a consistent background noise following -5/3 power law between the periods of 2 h and 4 days. Deseasonalized monthly data were used to analyse trends and extreme deviations (wavestorms). The possible trends, being obviously parts of longer variations, were not statistically significant over the studied 11 years.
The lists of most influential storms were different depending on the sites. The randomness of wave-storm events was higher in northerly exposed coasts (Letipea) where the longest fetches coincide with the northerly or northeasterly directions where strong storm winds are infrequent and quite random. On the other hand, westerly storm winds can be expected every winter on the west coast. A number of important European windstorms on course also reached Estonia and if their pressure minima made landfalls somewhere in southern Finland, prominent wavestorms were created in the Estonian coastal sea, too. The most influential wavestorm for the majority of the Estonian coast was the 9 January 2005 storm (Gudrun). A remarkably high (H s = 5.11 m) wavestorm was measured near Suurupi in December 2013.