Time-clustering of wave storms in the Mediterranean Sea

Giovanni Besio, Riccardo Briganti, Alessandro Romano, Lorenzo Mentaschi, and Paolo De Girolamo3 1Department of Civil, Chemical and Environmental Engineering, University of Genoa, Italy 2Department of Civil Engineering, University of Nottingham, UK 3Department of Civil, Architectural and Environmental Engineering, La Sapienza University, Rome, Italy 4European Commission, Joint Research Centre (JRC), Ispra, Italy Correspondence to: Giovanni Besio (giovanni.besio@unige.it)


Introduction
In recent years the occurrence of different coastal storms in a short time has been studied in the context of storm-driven erosion of beaches and dunes.Indeed it has been shown by different authors (Vousdoukas et al., 2012;Coco et al., 2014;Splinter et al., 2014;Karunarathna et al., 2014;Dissanayake et al., 2015) that storms occurring in quick succession may result in greater beach erosion than the cumulated erosion induced by single storms of far higher return periods.
In the events analysed in the aforementioned studies both the surge and the wave components played an important role.While studies that identify time clustering of storm surges are available (e.g.Wadey et al., 2014;Haigh et al., 2016), there is no study, to the best knowledge of the authors, that analyses the clustering properties of wave storms alone.In micro-tidal environments, such as the Mediterranean Sea, wave storms are the principal driver of short term coastal erosion and flooding; hence it is important to understand the occurrence of clustering.
The Mediterranean Sea wave climate has been extensively studied (e.g.Sartini et al., 2015a) and it is known that throughout the basin winter is richer in cyclones and, in turn, in wave storms.However, regional differences are significant.Sartini et al. (2015a) linked the seasonality of wave storms to local features of atmospheric pressure over the Mediterranean basin, strongly suggesting that the local typical meteorological conditions determine different temporal regimes of storm waves.
The present work addresses the gap in the knowledge of the occurrence of time clustering of wave storms by carrying out an analysis of wave storms sequences using the Allan factor (hereinafter AF, Allan, 1966;Barnes and Allan, 1966), a well-established technique to study the time behaviour of environmental processes.When the underlying process is characterized by clustering, the AF of a specific sequence of events is larger than 1 and shows a power-law behaviour at the timescales that exhibit departure from a Poisson distribution.The simplicity of the AF analysis made it popular in the study of time sequences of a number of physical processes Published by Copernicus Publications on behalf of the European Geosciences Union.such as earthquakes (Telesca et al., 2002;Cavers and Vasudevan, 2015), lightning (Telesca et al., 2008), rainfall (Telesca et al., 2007;García-Marín et al., 2008) or fires (Telesca and Pereira, 2010).However, the AF can also be larger than 1 for non-homogeneous Poisson processes, as shown in Serinaldi and Kilsby (2013).Hence it is important to distinguish clustering dynamics from cyclic Poisson processes.Methodologies that are suitable to achieve this are presented in Serinaldi and Kilsby (2013) and Telesca et al. (2012).
Here we analyse the AF on long time series of wave height in the Mediterranean Sea provided by hindcast reanalysis spanning the period 1979-2014 (Mentaschi et al., 2015).This analysis is validated and compared against the AF evaluated using the time series of wave measurements of the Italian national Sea Wave Measurement Network (Rete Ondametrica Nazionale, hereinafter RON).Subsequently we apply the methodology proposed in Serinaldi and Kilsby (2013) to gain an insight into the type of process that is described by the AF.The objective of this study is to identify the presence of time clustering of wave storms in the whole of the Mediterranean basin and examine the timescales at which events are correlated as well as the spatial distribution of the clustering.To this end, after scaling properties of wave storms are identified, they are mapped over the Mediterranean Sea.
The paper is organized as follows: after the Introduction, Sect. 2 explains the methodology used for the AF analysis, Sect. 3 describes the data sets used, Sect. 4 illustrates the results and Sect. 5 discusses the results and draws the conclusions of this work.

Clustering analysis methodology
Sequences of natural events such as earthquakes, rainfall and wildfires, can be seen as realizations of stochastic point processes.A process of this kind describes events that occur randomly in time and is completely defined by the times at which these events occur.Here time series of sea states are considered.Each sea state is defined by a set of spectral parameters, such as the significant wave height H s , the peak period T p , the mean period T m−1,0 and the mean direction of propagation θ m .Waves are always present on the sea surface; hence a sequence of storms needs to be extracted from a time series of sea states by only considering events that satisfy a certain criterion.A storm is commonly defined as a sequence of sea states in which H s exceeds a given threshold (e.g.Goda, 1988).In this work, a threshold for each node is defined by considering the local 98% percentile of the H s distribution, regardless of θ m (omnidirectional analysis; see Fig. 1 for threshold values of H s obtained with the hindcast model used here).The time t i at which the threshold is exceeded for the first time in each storm defines the event as part of a point process.If the interval between two subsequent events is below 12 h, the two are regarded as one event.This is common practice in analysing storms and the value is deemed appropriate for the Mediterranean Sea (e.g.Sartini et al., 2015a).Therefore, in each of the computational nodes over the Mediterranean Sea (see Fig. 2 for a map of the domain and the location of few control grid points used in this study to show the single-point behaviour of the AF), a point process is defined.An example for the control point A and for the years 2004 and 2005, is given in Fig. 3.In this figure it is evident that most of the storms during the 2 years considered occur between November and May, showing the pronounced seasonality that characterizes the basin.Figure 4 shows the number of events defined in each month over the year in the hindcast record for the same reference point A during the period 1979-2014 as a function of the percentile threshold (different wave heights).The seasonal variability of the storms in the Mediterranean basin is again recognizable.Note that the difference in the number of storms between the different percentiles considered is maximum in the most active months and, if the 99 % is chosen, the differences among seasons are small, although the seasonal variability is still recognizable.
These point processes are studied by defining equally spaced time windows of duration τ and counting the events in each window.The result is a sequence of counts N k (k = 1, . . ., M, where M is the number of time windows).The clustering of the events is then studied with the Allan factor (Allan, 1966;Barnes and Allan, 1966), defined as the variance of successive counts as In general terms, a point process is called fractal when a number of the relevant statistics shows scaling with related scaling exponents (Lowen and Teich, 1995).This implies that the AF depends on τ with a power law, with exponent α, which indicates the presence of clusters of points over a number of timescales τ .For a fractal process with 0 < α < 3 this power law (Telesca and Pereira, 2010) reads as follows: where τ 1 is the fractal onset time that marks the lower limit for significant scaling behaviour for the AF.For times smaller than τ 1 there is no significant time correlation, while for times greater than τ 1 a characteristic fractal trend can be derived from the value of the exponent.If the storm's process is Poissonian, the arrival times are uncorrelated; hence α is expected to be zero and the AF will be near unity.If non-Poissonian processes are present over a significant range of timescales it will be possible to identify α > 0 and AF > 1. Serinaldi and Kilsby (2013)   that cyclic, hence non-homogenous, Poisson processes show AF > 1 for timescales associated to cyclic components.It is therefore necessary to identify and separate the timescales at which clustering occurs from those at which the point process is Poissonian.To this end it is necessary to compare the AF pattern found in the wave time series with that of a process of known properties.A cyclic Poisson process is generated here with the same integrate and fire (IF) technique used in Serinaldi and Kilsby (2013).The cyclic components are selected by looking at the dominant harmonic components obtained with the Fourier analysis.The exponent α is estimated for the timescales at which the process is not Poissonian.Note that different ranges of τ can reveal different time scaling (clustering) of the same process through different slopes of Eq. ( 2) due to different kinds of forcing (Telesca and Pereira, 2010).
The occurrence of subsequent wave storms can be interpreted as a realization of stochastic temporal point process that could attain a clustered character when a number of its underlying features exhibit some scaling as a function of some scaling power law.The presence of such characteristics reveals that the process follows some kind of clustering in time (Lowen andTeich, 1995, 1996;Telesca et al., 2002;Telesca and Pereira, 2010).There are different statistical measures available in the literature to characterize the counting process of a general physical phenomena.In the present study we decided to employ the Allan factor (Allan, 1966;Lowen and Teich, 1996) thanks to the fact that it does not saturate the exponent α at unity as other indicators such as the Fano factor do.The Allan factor is defined as the variance of successive counts for a specific counting time window T divided by two times the mean number of counts in the same counting window For a fractal process the Allan factor recovers a power law of the type over an extended range of counting windows τ , where α is the so-called fractal exponent that for white noise time series attains values close to zero (i.e. the signal is characterized by the absence of time correlations, homogeneous Poissonian process), while for time-clustered processes it shows values greater than zero.τ 1 represents the fractal onset time and marks the lower limit for significant scaling behaviour for the Allan factor: for times smaller than τ 1 there is no significant time correlation, while for times greater than τ 1 a characteristic trend can be derived from the value of the exponent; furthermore different time windows can reveal different time scalings of the same process through different slopes of Eq. ( 4) due to different kind of forcing (Telesca and Pereira, 2010).Serinaldi and Kilsby (2013) demonstrated that cyclic, hence non-homogenous, Poisson processes show AF > 1 and power law behaviour for timescales associated to cyclic components.It is therefore necessary to compare the AF pattern found in the wave time series with that of a process of known properties.The same technique used in Serinaldi and Kilsby (2013) has been used here to simulate surrogate non-homogeneous point processes and compare them with the reference ones using a Monte Carlo approach.The processes have been generated using the same IF technique in Serinaldi and Kilsby (2013), to which the reader is referred for details.This further clarifies the nature of the process described by the AF and the role of the different cyclic components that contribute to generate above-threshold events.
3 Wave data

Wave hindcast
Wave hindcast in the Mediterranean Sea has been implemented on a time window covering 36 years, from the first of January 1979 to 31 December 2014 (http://www.dicca.unige.it/meteocean/hindcast.html).The wave model is forced by the 10 m wind fields obtained by means of the nonhydrostatic model WRF-ARW (Weather Research and Forecasting -Advanced Research WRF) version 3.3.1 (Skamarock et al., 2008).In the present study a Lambert conformal grid covering the whole of the Mediterranean Sea with a resolution of about 0.1 degree in longitude and latitude has been used.Initial and boundary conditions for atmospheric simulations were provided by the CFSR (Climate Forecast System Reanalysis) database (Saha et al., 2010).Use of CFSR reanalysis data for wave modelling provides reliable results, even if sometimes extreme wave conditions are not properly modelled (Cavaleri, 2009;Cox et al., 2011;Splinder et al., 2011;Carvalho et al., 2012;Chawla et al., 2013).For further details of the set-up and validation of the meteorological model, readers can refer to Cassola et al. (2015Cassola et al. ( , 2016)).
Generation and propagation of sea waves have been modelled using WavewatchIII ® , version 3.14 (Tolman, 2009).A 336 × 180 regular grid covers the whole of the Mediterranean Sea with a resolution of 0.1273 • × 0.09 • , corresponding to about 10 km at the latitude of 45 • N. Spectral resolution is characterized by 24 bins in direction and 25 frequencies ranging from 0.06 to 0.7 Hz with a step factor of 1.1.The output has been recorded hourly in all points of the computation grid for integrated quantities (i.e.significant wave height H s , mean period T m−1,0 , peak period T p , mean direction θ m , peak direction θ p , directional spreading θ ).The validation of the wave hindcast has been carried out through extensive comparison of simulated quantities and wave buoy data (see Mentaschi et al., 2013aMentaschi et al., , b, 2015) ) and has already been employed for different applications such as wave energy resource assessment (Besio et al., 2016) and extreme and wave climate analysis (Sartini et al., 2015a, b).

Buoy data
The Italian Sea Wave Measurement Network (Rete Ondametrica Nazionale RON) started operating in July 1989 (De Boni et al., 1992;Arena et al., 2001;Corsini et al., 2004).The locations of the buoys are indicated in Fig. 2. Until 1998 the network was made of eight pitch-roll directional buoys located offshore, in deep water conditions, of several sea areas equally spaced along the Italian peninsula.These original eight stations were La Spezia, Alghero, Ortona, Ponza, Monopoli, Crotone, Catania and Mazara del Vallo.The statistical wave parameters (i.e.significant wave height H s , mean period T m , peak period T p , mean direction θ m ) were originally retrieved every 3 h, below a station-dependent threshold for H s , and every half an hour above this threshold.The wave data time series, measured by the RON buoys, that have been analysed in the present study, cover a time window of 20 years, from the summer of 1989 to the spring of 2008 for the original eight buoys.For the cluster analysis performed using the RON records, data were considered every 3 h for all the stations.

Comparison between hindcast and buoy measurements
In order to assess the reliability of the hindcast time series related to storm cluster analysis, the results of AF for the RON buoys are analysed and compared to the corresponding grid points of the hindcast model.These results are shown in Figs.5-6.Results obtained on the basis of the RON data and hindcast series show a good qualitative and quantitative agreement, especially for lower threshold conditions (98 % percentile), while for higher thresholds (99.5 % percentile) they tend to present stronger differences, e.g. in Alghero (see Fig. 5).These findings can be explained by the fact that increasing the threshold limit would select just the most energetic wave conditions that are the most difficult to be reproduced by numerical models (a.o.Cavaleri, 2009) and sometimes to be recorded by wave buoys (breakdown, damages or even loss of the instrumentation).Also, differences are usually larger for smaller timescales, i.e. 0.5 < τ < 50 days and for the 99.5 % percentile (e.g.Alghero and Mazara in Fig. 5).These results confirm that the hindcast data and the wave buoys show very similar scaling properties.

Comparison with a simulated non-homogeneous point process
The AF patterns of both the model and data show a consistent behaviour across the Mediterranean basin.The AF is greater than one for τ greater than 12-24 h (0.5-1 days) and a distinct slope is recognizable, generally between 0.5 to 20-50 days at many of the stations.For larger values of τ , the AF increases to reach a maximum at 180 days.It is necessary to clarify the nature of the processes described by the AF patterns seen and, in particular, it is necessary to identify whether deviation from a cyclic Poisson process is present.
To this end, the AF pattern found from hindcast time series is compared with that of a simulated non-homogeneous Poisson process.This is generated using the IF technique em- ployed in Serinaldi and Kilsby (2013).The rate function of the simulated non-homogeneous Poisson process is generated as a sum of sinusoidal components with amplitudes, periods and phases obtained from the Fourier analysis of the reference signal.A Monte Carlo simulation of 1000 time series is then carried out and the simulated population of AF is compared with the reference one.Hindcast points A, G and O (see Fig. 2) are chosen for this analysis because they show different AF patterns in the timescales τ < 50 days.This analysis reveals that, as expected, the dominant cyclic component for all the considered time series is the one with a 1-year period.This was also noted for the RON data in Briganti and Beltrami (2008), where the amplitude of the annual cycle component was estimated to be around 0.25 m in Alghero, which is consistent with what was found in the present work.Together with the annual cycle the components with periods of 6, 3 and 1 months and 1 week have also been considered to simulate the non-homogeneous Poisson processes.The results of the comparison are shown in Fig. 7.For all three points it is clear that the simulated cyclic Poisson process explains the pattern of the AF at τ > 50 days well in all cases.As expected, this is the signature of the annual cycle, which strongly influences the occurrence of abovethreshold events.The AF departs from the Poisson distribution at τ < 50 days, above all in points A and G.The departure from Poissonian behaviour at these timescales occurs even at very low values of α, for example in point O.However, data often show oscillations, above all for α < 0.1, and it is not possible to make conclusions about the existence of a clustering regime.

AF results over the Mediterranean Sea
Results from the control points located over the basin (see    most cases.α varies significantly from point to point.A well-defined slope is very evident at points A (northern Tyrrhenian Sea), B (Gulf of Lion), D (Alboran Sea), and E (Algerian Sea).In all these cases a uniform value of α can be defined and the exponent value is in the interval 0.15-0.3.In other cases the slope is not so well defined or it is significantly smaller than 0.2.Points that show either or both characteristics are point R (Adriatic Sea), C (western Sardinia), F (Tunisian coast), G (southern Tyrrhenian Sea), M (Ionian Sea) and Q (Aegean Sea).At point Q (Aegean), α is virtually naught.
b.In the second group only the cyclic Poissonian regime is clearly recognizable, generally for τ > 20 days.At smaller scales the slope that is associated with the departure from the Poisson distribution is not present.This is the case for the southern Mediterranean points H (Egypt), I (western Libya), L (north-eastern Libya), O (south-eastern Mediterranean Sea) and P (southern Turkey).
The spatial distribution of the slope for small timescales is shown in Fig. 12.This figure has been obtained by determining the best fit value of α at different timescales.In order to take into account the local differences in determining the transition between slopes and the different regimes seen in the representative points, the slope has been estimated using four different ranges of τ .Clustering in the range 12 < τ < 72 h (3 days) is presented in Fig. 12a, for 12 < τ < 120 h (5 days) results are shown in Fig. 12b and finally Fig. 12c shows the results for 12 < τ < 240 h (10 days).Within this range the small-scale slope is higher in the northwestern Mediterranean Sea and, in particular in the northern Tyrrhenian Sea and in the Balearic Sea.Here α reaches values up to 0.3.Areas with α around 0.2 are present in the Adriatic Sea, on the Syrian and Lebanese coast and along the Tunisian coast.The effect of widening the range of τ is to decrease the best fit value of α.This effect reduces the regions that show α significantly higher than zero, in particular in the Adriatic Sea and on the eastern coast of Tunisia.When the interval 12 < τ < 240 h (0.5-10 days) is used (Fig. 12c), the best fit of α is significantly higher than zero only in the northwestern Mediterranean Sea with the average α around 0.2 and zones with α > 0 are present in the eastern part of the Adriatic Sea and on the Syrian coast.

Discussion and conclusions
The results presented highlighted the presence of a departure from the Poisson distribution for timescales shorter than τ < 1200 h (50 days).This regime is characterized by α = 0.15-0.3and is more evident in the north-west of the Mediterranean Sea.In the rest of the basin α is closer to zero and the AF pattern is characterized by oscillations, without a well defined regime.For τ > 50 days the arrival of above-threshold storms is dominated by the effect of seasonal and interseasonal oscillations and can be described as a cyclic Poisson process.Similar scaling regimes have been observed in other phenomena with seasonal behaviour, e.g.fires (Telesca and Pereira, 2010).These results match with the findings by Sartini et al. (2015a), who found that the northern basin RON buoys (e.g.Ponza and La Spezia buoys in the Tyrrhenian Sea) showed lower seasonality than the buoys in the southern basin (e.g.Crotone, in the Ionian Sea).La Spezia buoy, for example, is located in the Ligurian Sea, a region where departure from the Poisson distribution is higher.Although in the region the cyclogenesis in the Gulf of Genoa shows marked seasonality, cyclones are present throughout the year (Lionello et al., 2006;Sartini et al., 2015a).This persistence of cyclonic events helps to explain the behaviour at smaller scales (i.e.τ < 1200 h, 50 days).The clustering at scales of days indicates that meteorological conditions favour the occurrence of multiple events over a few days.It is not the case that this behaviour is seen in the most active cyclonic region of the Mediterranean Sea, e.g. the north-west according to Lionello et al. (2016).Similar considerations apply to the northern Adriatic Sea.In other parts of the basin, where these persistent conditions do not occur, the arrival of storms is well described as a cyclic-Poisson process.
The values of α found in the present study do not allow us to draw conclusions on whether this deviation from a Poisson distribution is large or small for the phenomenon at hand, as there is no comparison with other basins.Because of this, it is important to analyse other basins.
The clustering at the timescales found has the potential to exacerbate local beach erosion generated by individual storms, as shown in Dissanayake et al. (2015); hence it will be important to understand the implication of these time regimes on the dynamics of the Mediterranean coastal regions.

G
Figure 1.Value of significant wave height threshold in metres for the 98 % percentile.

Figure 2 .
Figure 2. Hindcast control grid points (red circle) and RON buoys as reference points (yellow circles).

Figure 3 .
Figure 3. Storm occurrence for the northern Tyrrhenian reference point (A): 2004/2005 in the top panel, zoomed-in graph of winter 2004/2005 in the bottom panel.

Figure 4 .
Figure 4. Number of storms vs. threshold for the northern Tyrrhenian reference point (A).

Figure 5 .
Figure 5.Comparison of Allan factor between RON and hindcast data series for different threshold percentiles (98 and 99.5 %).

Figure 6 .
Figure 6.Comparison of Allan factor between RON and hindcast data series for different threshold percentiles (98 and 99.5 %).

Figure 7 .
Figure 7.Comparison of Allan factor between hindcast data series for 98 % percentile (black line) and 1000 simulated cyclic Poisson processes (grey lines).The AF corresponding to the 95 % percentile of the AF distribution is also plotted (dashed line).Top left shows point A (northern Tyrrhenian).Top right shows point G (southern Tyrrhenian).Bottom shows point O (south-eastern Mediterranean).

Fig. 2 )
are shown in Figs.8-11.The analysis of the AF curves reveal that these can be divided in two groups: a.The first group clearly shows the slope corresponding to the departure from the Poisson regimes.The change in regimes occurs at around τ = 50 days in

Figure 8 .
Figure 8. Allan factor (AF) as a function of counting window τ and of the wave height threshold (different percentiles as in the legend) for different locations in the Mediterranean Sea (cf.Fig. 2).

Figure 9 .
Figure 9. Allan factor (AF) as a function of counting window τ and of the wave height threshold (different percentiles as in the legend) for different locations in the Mediterranean Sea (cf.Fig. 2).

Figure 10 .
Figure 10.Allan factor (AF) as a function of counting window τ and of the wave height threshold (different percentiles as in the legend) for different locations in the Mediterranean Sea (cf.Fig. 2).

Figure 11 .
Figure 11.Allan factor (AF) as a function of counting window τ and of the wave height threshold (different percentiles as in the legend) for different locations in the Mediterranean Sea (cf.Fig. 2).

Figure 12 .
Figure 12.Spatial distribution of the exponent α for the whole of the Mediterranean basin. demonstrated