**1. Introduction**

Coastal areas are extremely important in Mediterranean countries since they host the majority of their population and economic activities [1]. Over the last few decades, one of the faster urban developments has occurred along the Spanish Mediterranean coast, especially at the Costa del Sol [2]. As a result of this expansion, human activities and buildings were placed extremely close to the shore [3]; therefore, they are now threatened by natural hazards influenced by climate change-related processes such as sea-level rise and increases in storm frequency and intensity [4,5]. To reduce storm impacts, it is necessary to understand specific coastal characteristics and sensibilities as well as to fully comprehend storm nature. In recent years, several researchers have studied these aspects from different viewpoints. In Spain, Rodrıguez-Ramırez et al. [6] studied storm records on the Huelva coast to obtain appropriate future development and management strategies; Mendoza and Jiménez [7] and Mendoza et al. [8] presented an intensity scale for wave storms on the Catalan coast to characterize their spatial and temporal variability; Guisado and Malvárez [9] and Pintado and García [10] used extreme wave conditions to complete the characterization of the morpho-dynamic environments of the Costa del Sol and Huelva areas; Anfuso et al. [11] and [12,13] characterized storms along the Atlantic side of Andalusia. In recent decades, coastal scientists have used several indexes to characterize storms, e.g., Halsey [14] ranked north-east Atlantic coastal storms (northeasters or nor'easters) based on a damage potential index and Dolan and Davis [15] proposed an intensity scale index to classify nor'easters into five classes, from weak to extreme, based on wave height and storm duration. Orford et al. [16] and Orford and Carter [17] used the role of storm surge to develop a new storm index. Kriebel et al. [18] proposed a nor'easter risk index by combining the effects of storm surge, wave, and duration and

Zhang et al. [19] developed a storm erosion potential index by combining the effect of storm tide, wave energy, and duration. This paper analyzes a 35-year wave climate dataset obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) for four available prediction points equally spaced along the Mediterranean coast of Andalusia (south of Spain). This allowed the definition and assessment of storm characteristics and their spatial and temporal distribution along the investigated area. To characterize the storms, a new approach was adopted, assessing the real wave energy flux of each storm, using a robust definition of the storm itself. During the investigated period, a total of 2961 storm events were recorded. These were classified according to five classes of storms, from low (Class I) to high-energetic (Class V). Results obtained are useful to understand potential impacts of both single and grouped storms, and hence put in place the appropriate prevention and mitigation strategies.

### **2. The Study Area**

This study is focused on the wave climate of the Andalusia Mediterranean coast, a very populated area whose land cover has experienced important changes during recent decades [20]. Málaga is the province that has experienced the most important coastal occupation, in particular due to the construction of structures related to national and international tourism [20]. The coast is a micro-tidal environment (tidal range < 20 cm, [9] ), about 546 km long, and including four provinces, i.e., Cádiz, Málaga, Granada, and Almería (Figure 1).

**Figure 1.** Geographic location of the study area and wave rise (1979–2014) for each prediction point.

There are many Andalusian coastal areas that suffer erosion problems [21], typically linked to very energetic storms producing severe damage to coastal structures. For example, Figure 2 shows the damage produced in winter 2015 by a western storm at Alumuñecar (Granada Province). This storm particularly affected the beaches of San Cristobal and La Herradura. At the former, the extreme wave run-up broke the facilities for summer tourism, and at the latter, storm waves reached the road at several points, depositing cobbles and sand that endangered people the circulation of vehicles.

**Figure 2.** Severe damage after a western storm in winter 2015 at the Reina Sofia Promenade (Almuñecar), photo by europapress.

On the Mediterranean Andalusian coast, the beaches are rectilinear and composed of medium to fine and dark to golden sands and, especially in Granada and Almería, cliffed sectors are observed. The near-shore area generally shows high slope values and intermediate and reflective morpho-dynamic states linked to a narrow continental shelf [22]. This sector is exposed to winds from E and SE with minimum and maximum wind speed values that range from 0.4 to 0.9 m/s [9].

#### **3. Methods**

With the aim of characterizing the wave climate of the studied area, four prediction points were identified along the Mediterranean Andalusian coast; at each point, wave rises and the monthly means of maximum wave height *Hm*0,*max* were calculated. The wave rises give information about the direction and intensity of incoming waves, whereas the monthly wave height means give information about the seasonal characteristics of the Mediterranean Andalusian coast. To identify each single storm, the definition of [23] was adopted, which allowed calculation of the energy flux by using the linear deep-water wave theory and, finally, classification of energy flux, preferring this parameter to empirical ones.

## *3.1. Wave Climate Preliminary Analysis*

Wave climate analysis was carried out using wave data modelled by the ECMWF by means of the WAve Model (WAM). This numerical model, which solves the energy balance equation, forecasts wave climate that is then subjected to quality controls ensuring its consistency [24] (https://www.ecmwf. int/en/elibrary/16951-wave-model accessed on January 2019). This paper used MATLAB scripts to analyze wave characteristics along the Mediterranean coast of Andalusia for the period 1979–2014, obtained by the ECMWF within the framework of the ERA-INTERIM project. The four predictions points, from W to E, used in this paper, are Point 1, close the Strait of Gibraltar, Point 2, east of Málaga, and Points 3 and 4 in front and east of Almería, respectively (Figure 1). Each temporal series is formed by 51,860 data points recorded over a period of 35 years (1 January 1979–30 January 2014) which is enough to analyze potential trends of increasing wave heights, the presence of climate-controlled cycles, or annual variations due to climate events [5].

#### *3.2. Storm Classification Using the Energy Flux*

To determine storm events, the definition used was the one given by [23], i.e., a Mediterranean Sea storm is a sequence of sea states in which the spectrally significant wave height exceeds the threshold *ht* and does not fall below this threshold for a continuous time interval greater than 12 h. He also considered that the time interval between consecutive single storms must be greater than 12 h. The government agency "Puertos del Estado" [25] suggested, for the Mediterranean coast of Andalusia, adopting a threshold value *ht* = 1.5 m. Boccotti [23] suggests adopting the same value for the threshold (*ht*) because it is 1.5 times the mean yearly significant wave height. Following the aforementioned criteria, 2961 storms were selected for the period 1979–2014.

Many studies (e.g., [7,8,11–13]) based the classification of storms on the use of the Dolan and Davis [15] "Storm Power Index". In this paper, a physically based parameter was preferred, namely, wave energy flux [26,27]. Wave energy flux, or wave power per unit of wave-front length (*P*), was calculated using the following equation:

$$P = \frac{\rho g^2}{64\pi} T\_\epsilon H\_{m\_0}^2 \left[\frac{W}{m}\right] \tag{1}$$

where *ρ* is water density, *g* is the gravity acceleration, *Te* is the energy period that represents the period of the sinusoidal wave with the same energy as a real sea state (for which a JONSWAP spectrum is about 90% of the peak period) and *Hm*<sup>0</sup> is the spectrally significant wave height. To obtain an accurate estimation of the total energy (*E<sup>i</sup> tot*) of each storm [26–29] the energy flux was time-integrated:

$$E\_{tot}^i = \int\_0^{d\_i} Pdt \left[\frac{Wh}{m}\right] \tag{2}$$

where the *d<sup>i</sup>* is the duration of *i*-th storm. Using Equations (1) and (2) the total energy of each of the 2961 storms was calculated.

For example, two storms at the prediction Point P3 are shown in Figure 3. The first storm started on 25 February 2009 and the second on 4 March 2009. Figure 3a shows the *Hm*<sup>0</sup> values of the two storms identified by means of the threshold *ht* = 1.5 m. The corresponding energy flux *P* and the total energy *E<sup>i</sup> tot* (Equations (1) and (2)) is shown in Figure 3b.

**Figure 3.** (**a**) Spectrally significant wave height *Hm*0, during the storms of 25 February and 4 March 2009 at the prediction Point P3. The dashed red line is the spectrally significant wave height threshold *ht* = 1.5 m, *d*<sup>1</sup> and *d*<sup>2</sup> are storm duration. (**b**) Energy flux during time and total storm energy.

The method proposed by [30], i.e., "the natural breaks" function, was used to classify storm events into five classes, from Class I (low-energetic events) to Class V (high-energetic events).

#### The Return Period of the Energy Flux of Storms

For the estimation of the return period of the energy flux (*P*) of each storm class, the probability of exceedance was fitted with a modification of the Weibull Cumulative Distribution Function (CDF) described in [23,31]:

$$F(P) = e^{-\left(\frac{p}{w}\right)^{\mu}} \left[\%\right] \tag{3}$$

where *w* and *u* parameters that depend on the location under examination. Using the auxiliary variables *X* = 100 · *ln*(2.5 · *P*) and *Y* = 100 · *ln* [*ln* (1/*F*(*P*))] [23,31], in the coordinate system *X*, *Y*, the data points should lie on a straight line. Thus, the parameters of Equation (3) can be easily estimated by fitting those data points by means of the least-squares method.

The above-mentioned Weibull CDF is related to the return period *Tr* by means of

$$T\_r = \frac{1}{\lambda \cdot F(P)} \left[ year \right] \tag{4}$$

where *λ* is the mean number of events per year. Consequently, for the lower and upper limit of each storm class, the probability of exceedance and the return period were calculated.

#### *3.3. Stormy Year*

For the sake of continuity with previous research [11,32], the definition of stormy year was adopted here. In particular, for each prediction point, stormy year empirical recurrence period (*Ri*) and annual frequency of occurrence (*fo*) were assessed by using the equations:

$$R\_i = \frac{(n+1)}{m} \left[ year \right] \tag{5}$$

$$f\_0 = \frac{1}{R\_i} \begin{bmatrix} \% \end{bmatrix} \tag{6}$$

where *Ri* is the recurrence interval calculated for *n* number of years (35 in this paper), *m* is the number of events that occurred within the date-range of interest, and *fo* is the yearly frequency of occurrence of the event. For each prediction point, two values of *Ri* and *fo* were calculated using a minimum and a maximum threshold, respectively. The minimum threshold is the mean of total energy (*μ*) calculated within the whole period of 35 years, the maximum threshold is *μ* + *σ*, where *σ* is the standard deviation of the total energy.

#### **4. Results and Discussion**

#### *4.1. Wave Height Characterization*

Wave height data did not present a general trend along the investigated period, but a clear seasonal behavior was recognized, as observed by Rangel-Buitrago and Anfuso [13] on the Atlantic sector of Andalusia. As a general trend, higher average values of monthly maximum spectrally significant wave heights (*Hm*0,*max*) were observed during the winter season, i.e., the December–March period and, specifically, Points 1, 2 and 4 recorded maximum values in March and Point 3 in February (Figure 4). High values recorded in March are linked to the great importance of eastern waves due to regnant winds during such months. Lower average *Hm*0,*max* values were observed during the summer time, i.e., July, August, and September, ranging from 1.7 m at Point 1 to 2.1 m at Point 3. Dealing with the behavior of each analyzed point, Point 3 showed the highest values in all years, always followed by Point 2; but during June, July, and August, Point 4 recorded greater values than Point 2 (Figure 4).

**Figure 4.** Average values of monthly maximum spectral significant wave heights *Hm*0,*max*.

Approaching directions observed at different points clearly reflected coastal orientation and prevailing marine climate (Figure 1). Point 1 is close to the Strait of Gibraltar but sheltered to the Atlantic swell waves. East approaching fronts, with prevailing wave height classes of 0.5–1.25 m (c. 30%) and 1.25–2.5 (c. 10%), clearly prevailed at this prediction point (Figure 1) and were linked to the strong easterly winds associated with a surface pressure gradient over the Gibraltar Strait when the Azores high pressure is located over the Iberian Peninsula, while there is pronounced low pressure over northern Africa [33]. Points 2 and 3 are situated at the central part of the investigated area, so they are exposed to winds and waves from W and E and E-NE directions; specifically, at Point 2, E and N approaching waves are equivalent to W approaching fronts. At Point 3, the western component prevailed on the north-east approaching direction because of the increase of the western geographic fetch; an increment of the NE component was also observed since this point is more exposed to this approaching direction with respect to Point 2 because of coastal orientation (Figure 1). At Point 4, despite the prevalence of the E-NE approaching direction, the NE component becomes even more evident than at Point 3. Furthermore, since this latter point is sheltered to the west because of coastal orientation (NNE-SSW oriented), waves approaching from the third quadrant essentially present SW and W-SW components.

#### *4.2. Storm Characterization*

During the investigated period, a total of 2961 storm events were categorized into five classes, i.e., Class I (weak), Class II (moderate), Class III (significant), Class IV (severe), and Class V (extreme). Points 3 and 2 recorded the highest number of storms (Figure 5). The distribution is similar at all points, with a clear dominance of events belonging to Classes I and II: 87.4% in Point 1, 86.6% in Point 2, 83.1% in Point 3, and 90.4% in Point 4 (Table 1, Figure 5). Mean wave height value of each class did not present great spatial variations, i.e., all points recorded similar values of wave height for the same class. Concerning maximum and minimum wave height values per class, a clear and general trend was not observed, e.g., Point 4 presented the lowest wave height value for Class III but the highest for Classes IV and V (see Table 1).


**Table 1.** Storm characteristics at each prediction point: class, energy (*Etot*), frequency of storms, significant wave height (*Hm*0) peak period (*Tp*), and duration (*D*).

**Figure 5.** Distribution of storm events per class at each point.

Regarding wave period, the same storm classes of different points presented similar values. Higher wave period variations among classes were recorded at Points 2 and 3 (Table 1). The mean storm duration presented an increase from Class I to Class V at each point (especially Points 2 and 3) and important differences among points (Table 1). Figure 6 reports the monthly distribution of all storms (i.e., the sum of all events recorded at each points) per class. Results are very similar to those obtained by [7,8,11–13,15,34] in their respective studies. Concerning temporal distribution, Classes I and II storms were observed along the whole year. Class III storms were recorded in winter and spring seasons (from October to May), with a minimal occurrence during summer months, i.e., June (8 storms), August (2 storms), and September (2 storms). Class IV storms were observed from November to March, and Class V storms were only recorded from December to March (especially in February, with 7, and March, with 5 events). It is interesting to observe that Classes II and III are very frequent in March, because it relates to approaching waves generated by E and SE winds that are quite frequent in springtime.

**Figure 6.** Monthly distribution of all storm events per class and month.

Concerning approaching directions, at Point 1, Classes I and II mainly approached from the east with a small western component (Figure 7) and Classes III, IV and V almost exclusively (>95% of records for each class) approached from E and were linked to the predominant eastern storm waves (Figure 1). At Point 2, the sum of the E and E-NE directions is broadly equivalent to the western component (Figure 7) and clearly reflects, with a slight increase of the eastern component, the wave rise shown in Figure 1. At Point 3, storm energy classes approaching directions (Figure 7) broadly reflect wave rise presented in Figure 1 with an increase of the importance of the western component for all storm classes. At Point 4 (Figure 7), the storm approaching directions are consistent with the wave direction pattern observed in Figure 1. At Point 4, more numerous storms approach from the E-NE direction and disappear in the SW direction. The return period (*Tr*) of Classes III, IV, and V events at each prediction point are shown in Table 2.

**Figure 7.** Energy flux roses at all prediction points. (**A**) Point P1, (**B**) Point P2, (**C**) Point P3, (**D**) Point P4. *N* is the total number of storm events.

At Points 2 and 3 the yearly the probability of energy flux exceedance for Class III storms was 100%, varying from 76.9% to 100% at Point 1 and from 43.5% to 100% at Point 4 (Table 2). For Class IV, the probability of occurrence ranged from 41.7% to 100% and from 58.8% to 100% at Points 2 and 3 respectively, and from 27.8% to 76.9% at Point 1 and from 12.7% to 43.5% at Point 4. Class V probability of exceedance have minimum percentages of 12.7% and 18.5% at Points 2 and 3, and 8.4% at Point 1 and 3.3% at Point 4. Mentioned values are broadly similar to observations carried out by Anfuso et al. [11] near the areas of Huelva and Cádiz, i.e., the less energetic zones of Cádiz Gulf. It is observed that the most energetic points (i.e., Points 2 and 3) have the highest percentages of probability of energy flux exceedance (Figure 8), so the facing littoral, i.e., the coast between Málaga and Almería, is very exposed to storms belonging to all classes, and especially to most energetic ones (III, IV, and V) that can have a great impact on both natural and urbanized sectors.

**Figure 8.** Energy flux probability of exceedance (%) of storms per class III, IV, and V at each point in the coastal area of influence. The percentage values per class as reported on Table 2 are superimposed in the colored stripes.

### *4.3. Characterization of Stormy Years*

In this paper, a special attention was devoted to the yearly cumulative characteristics of storms, i.e., their number, duration, and energy. In fact, as observed by Ferreira [35], the relationship between storms and beach erosion (and/or damage to human structures) varies according to single storm characteristics, storm grouping, and coastal response/morpho-dynamic behavior. Concerning the characterization/quantification of stormy years, only Classes III, IV, and V were considered since these events are the ones that produce important coastal damage according to [12]. With respect to energy data, a similar trend was recorded at the four points, with a similarity observed among Points 1, 2, and 3 (Figure 9).

**Figure 9.** In the left column: the total energy of storms per year for Classes III, IV, and V. At (**a**) Point P1, (**c**) Point P2, (**e**) Point P3, and (**g**) Point P4. In the right column: the corresponding durations of storms (**b**) Point P1, (**d**) Point P2, (**f**) Point P3, and (**h**) Point P4. Red and dashed lines respectively represent the average value (*μ*) and the average plus one standard deviation (*μ* + *σ*).

Dealing with yearly energy distribution, it was possible to observe nine energetic years, i.e., 1980, 1983, 1990, 1992, 1995, 2001, 2008, 2010, and 2013. Anomalies were recorded at Point 3, which differed in four of the nine aforementioned years; at Point 1, which recorded a low-energy year in 2013 (i.e., essentially Classes I, II, and III), and at Point 4, which recorded low energy values in 1983, 1984, and 2008 and high values in 1986 and 2007 (Figure 9). Analyzed data of stormy years did not present a clear trend; as an example, yearly distribution of cumulative energy presented a correlation factor (*r*2) that ranged from 10−<sup>5</sup> for Point P2 to 2 × <sup>10</sup>−<sup>4</sup> for Point P4, but showed a cyclical behavior as highlighted by the calculation of the return period. The recurrence interval and the yearly frequency of occurrence are presented in Table 3 according to the distribution of yearly cumulative energy values (Figure 9). Per each Point, two values of recurrence interval (and frequency of occurrence) are presented, and refer to the mean value and the mean plus one standard deviation of yearly cumulative energy. The frequency of occurrence of stormy year was higher at Point 2 with 17% and 53%, respectively, for the lower and higher energy value, i.e., a low-energetic storm year is observed every 1.9 years and a high-energetic year presents a recurrence interval of 6 years. Point 4 recorded longer recurrence intervals, which ranged from 2.6 to 9 years (39% and 11% values of frequency), respectively, for higher and lower energetic stormy years. Points 1 and 3 presented average values (Table 3).


**Table 2.** Return period (*Tr*) and probability of energy flux exceedance (*F*(*P*)) for the lower and upper limit of each Class (III, IV, V) calculated using Weibull CDF Equation (3). The minimum values are highlighted with blue-colored font.

**Table 3.** Stormy year recurrence interval (*Ri*) and yearly frequency of occurrence (*fo*) calculated from Equation (4), considering yearly cumulative energy values (Figure 9), for each prediction point (Table 1).


Regarding the cumulative storm duration (Figure 9), the trend was not clear. From one side, there is a good correspondence between the number of storms and storm duration within each point, that is, years with many storms also recorded a high cumulative yearly storm duration. Such correspondence was also observed with yearly energy values distribution. On the other hand, no correspondence was observed among different points. A comparison of stormy year distribution recorded in this paper with the one observed by [12] for the Atlantic region of Andalusia showed the opposite behavior linked to the fact that storms on the Atlantic side of Andalusia are essentially related to the predominance of western approaching fronts; meanwhile, on the Mediterranean side [36], they are essentially linked to eastern approaching fronts.
