Trends and variability of ocean waves under RCP8.5 emission scenario in the Mediterranean Sea

Wind-generated ocean waves are key inputs for several studies and applications, both near the coast (coastal vulnerability assessment, coastal structures design, harbor operativity) and off-shore (a.o. oil and gas production, ship routes, and navigation safety). As such, the evaluation of trends in future wave climate is fundamental for the development of efficient policies in the framework of climate change adaptation and mitigation measures. This study focuses on the Mediterranean Sea, an area of primary interest, since it plays a crucial role in the worldwide maritime transport and it is highly populated along all its coasts. We perform an analysis of wave climate changes using an ensemble of 7 models under emission scenario RCP8.5, over the entire Mediterranean basin. Future projections of wave climate and their variability are analyzed taking into account annual statistics of wave parameters, such as significant wave height, mean period, and mean direction. The results show, on average, a decreasing trend of significant wave height and mean period, while the wave directions may be characterized by a slight eastward shift.

busiest seas in the world (Leone 2017). Therefore, a detailed knowledge and monitoring of oceanographic variables in this basin is of foremost importance. In particular, gravity waves play a relevant role on several aspects of ocean and coastal dynamics. They are one of the main drivers of coastal erosion and accretion (De Leo et al. 2016;Harley et al. 2017;Mentaschi et al. 2018). Extreme and multimodal sea states can pose a threat to the safety of navigation and to coastal structures (Soares and Teixeira 2001;de Osés and Castells 2008;Ventikos et al. 2018). Wave setup and runup contribute to the extreme sea levels that drive coastal floods (Didier et al. 2015;De Leo et al. 2019) and coastal erosion (Shih et al. 1995;Ruggiero et al. 2001), together with storm surges and tides (e.g., Cazenave and Cozannet 2014;Melet et al. 2018).
According to projections of future climate, the ongoing sea level rise (SLR) will increase the magnitude and frequency of coastal hazard worldwide (Hinkel et al. 2014;Vousdoukas et al. 2018). However, local increases/decreases in wave height and periods or directional shifts of waves could exacerbate or mitigate this tendency. The important question of understanding how the wave climate will evolve in view of future climate changes has been tackled by a number of research groups worldwide, which in the last decade started coordinating their efforts within the Coordinated Ocean Waves Climate Project (COWCLIP, www.cowclip.org). Among the existing global and largescale studies, it is worth mentioning the joint COWCLIP analyses Morim et al. 2019;Morim et al. 2020a), as well as the publications of individual research groups on mean wave climate and on the extremes (e.g., Mori et al. 2010;Semedo et al. 2012;Wang et al. 2014;Perez et al. 2015;Shimura et al. 2016;Aarnes et al. 2017;Mentaschi et al. 2017;Bricheno and Wolf 2018). On the basis of global altimetry Young and Ribal (2019), and earlier Young et al. (2011), produced analyses of the global historical trends of significant wave height (H s ). Other authors estimated how changes in wave climate will translate into wave setup and wave runup along the world's coasts Melet et al. 2019). These studies shed light on the present and future changes of wave climate and on their consequences; however, they lack the spatial resolution for a comprehensive analysis in smaller basins or at a regional scale. Few studies were developed specifically for the MS or its sub-basins, notably the work of De Leo et al. (2020), the pioneering work by Lionello et al. (2008), based on a single model with a resolution of 50 km, and the study by Casas-Prat and Sierra (2013), limited to the northwestern Mediterranean. The advancement of our knowledge of local basins is the objective of a specific task of the COWCLIP initiative, dedicated to the development and intercomparison of comprehensive regional projections.
This contribution aims to improve our understanding of the future changes of wave climate in the MS, within the framework of the regional projection task of the COWCLIP initiative. We carried out wave simulations with the model Wavewatch III (WW3DG 2019) with a resolution of about 10 km in longitude and latitude, forced by an ensemble of 7 Euro-CORDEX regional models (Jacob et al. 2014) in Representative Concentration Pathway (RCP) scenario RCP8.5 (Van Vuuren et al. 2011). Then, we computed some annual statistics of the expected future significant wave height {H s }, mean period {T m }, and mean direction {θ m }, and studied the relative trends. Section 2 describes the regional climate models (RCM) and wave model used to derive the projections of future wave climate in the Mediterranean, along with the analysis employed to characterize trends in wave parameters. Results are outlined and discussed in Section 3, followed by conclusions in Section 4.

CORDEX-Forced regional wave model
Simulations for future wave climate conditions in the MS were performed employing the third-generation wave model WavewatchIII v5.16 with a resolution of 0.127 • in longitude and 0.09 • in latitude, corresponding to about 10 km, and a time output for all the wave bulk variables equal to 3 h. The model was implemented with the ST4 source term (Ardhuin et al. 2010), and a setup similar to that employed by Mentaschi et al. (2015) for the implementation and validation of a wave hindcast in the MS.
The future projections of wave climate were forced by seven different Euro-CORDEX climate projections under scenario RCP8.5 (see Table 1). These RCM were chosen due to their relatively high resolution (11 km), which is needed to resolve adequately the complex orography and geography of the MS. A similar assumption motivated the studies of Bricheno et al. (2013) and Menendez et al. (2014), who showed that an improved accuracy is achieved in coastal waters and in the MS owing to the regional downscaling of wind forcing. Due to the uncertainties inherent in modeling future climate, we opted to develop an ensemble approach employing seven different RCM in order to exhaustively investigate the variability of future projections of wave climate. Simulations available for control (historical) period and RCP8.5 scenario were carried out for the time slices 1970-2005 and 2006-2100, respectively, saving the maps (i.e., gridded fields) of H s , T m , and θ m every 3 h along with other variables over the entire MS. The historical output of each ensemble member, for each of the analyzed variables, was compared with hindcast values (Mentaschi et al. 2013(Mentaschi et al. , 2015 to assess the performance of the models. The hindcast was developed on the same computational grid used for the RCM projections, and provides hourly data of wave parameters in the 1979-2019 period. The common time slice 1979-2005 was used for the comparison: first, the time series of H s and T m were averaged over the whole period, resulting in a single data point for each node. Then, for each ensemble member, we computed the maps of the bias and the spatial correlation among the nodes between the hindcast and the RCM data.

Trend analysis
The annual means, annual 90th percentiles, and annual maxima of H s and T m were computed for each grid point using the COWCLIP utility getStat. The annual means of θ m were computed instead with the utility getStatDir (Morim et al. 2020b , respectively. A similar notation is used for the time series of T m statistics. The yearly ensemble values were then computed as the mean over the seven models. As for θ m , the arithmetic mean was first computed on the components projected into Cartesian coordinates, then it was subsequently converted back to polar coordinates  (Fisher 1995). Such an approach allows to account for the discontinuity in the 0/2π space: given a series of n directions θ m j , the circular mean can be defined as: In case of the θ m time series, the computation of the ensenmble mean values took advantage of the CircStat Matlab toolbox (Berens et al. 2009).
Then, analyses of wave climate changes were performed on the 2006-2100 time series, that is, for the future time slice. First, trends on the time series of H s and T m annual statistics were quantified through the slope of the respective best linear fit. In this work, we employed the Theil-Sen's slope (Sen 1968, Theil 1992, henceforth referred to as b. Given a series of x i data (i = 1...n, n being the number of samples), its computation reads: where x j and x l are the j th and l th data of the series, respectively. Then, the (1-α) confidence interval of b was computed according to Hollander et al. (2013): where denotes the inverse of the standard normal distribution, and N is the number of slopes related to the n(n − 1)/2 possible pairs of sample points; this computation was performed through the Matlab builtin function provided by Burkey (2006). The confidence interval of b allows assessing whether the computed trend is significant or not. Namely, if a change of sign occurs between b + (the upper confidence level) and b − (the lower confidence level), this may indicate that the data are too dispersed to compute a reliable linear trend.
As a second step, the estimates of b were analyzed in the context of the Mann-Kendall test (Mann 1945;Kendall 1955). This test allows to check whether an either upward or downward monotonic trend is present within a dataset on the basis of the test p value, computed as: said p MK the p value related to the Mann-Kendall test, and Z MK being the test statistic, that depends on the ranks of the data within the series they belong to. The values of p MK can be used to evaluate whether the data behave consistently with the hypothesis that there is no trend (p MK close to 1) or, at the opposite, that a trend affects the data (p MK close to 0). Such a combined use of b and p MK was already used by De Leo et al. (2020) that showed how to use these information jointly to strengthen the analysis of long-term trends on time series of H s data. Furthermore, the present work referred to the innovative trend analysis (referred to as ITA, Ş en 2011, 2013), a simple and intuitive graphical analysis. Given a time series of data, the ITA implies selecting two subsets of equal length belonging to successive time periods. These subsets are subsequently sorted in ascending order, and plotted versus each other in a common range. If the resulting scatters lie entirely above (below) the 45 • straight line, the existence of a positive (negative) trend can be deduced. An example  of the ITA for three synthetic datasets characterized by different trends is shown in Fig. 1. In this case, the sorted subsets are referred to as x 1 (belonging to the earlier time slice) and x 2 (belonging to the later time slice).
Here, the annual statistics of H s and T m were averaged across all the computational nodes, leading to the regional mean annual statistics of the parameters at hand. Then, the data of 2010-2040 and 2070-2100 time slices were selected and compared according to the ITA, to get further indications on the future trends of H s and T m . Similarly to De Leo et al. (2020), the trend magnitude resulting from the ITA was assessed computing the distances (hereinafter referred to as δ) of the scatters with respect to the 45 • straight line (which identifies the absence of trend). Positive and negative trends can be therefore discussed in terms of the empirical cumulative distribution function (ecdf) of delta, and of its position relative to the 45 • line. The ITA should be carefully used, as to rely only on this test could lead to wrong conclusions (Serinaldi et al. 2020). However, in this work, the outcomes related to three different tests were compared, in order to discuss and strengthen the outcomes related to each single method.
Besides the statistical significance of the trend, in a study on climate projections it is important to assess the consistency of the projected changes across the ensemble (Knutti and Sedláček 2013). In this work, this has been accomplished by evaluating the number of ensemble members indicating either positive or negative trends at each point of the computational grid. Finally, θ mean m was investigated by means of polar plots of the time series at six locations, selected in different sub-basins of the MS.

Analysis of historical wave climate simulations
First, the historical wave climate in the MS driven by each of the RCM was compared to the hindcast data developed and validated by Mentaschi et al. (2013) and Mentaschi et al. (2015). Attached in the Appendix are the comparison maps for H s (Figs. 19,20,and 21) and T m (Figs. 22,23,and 24), which show the different spatial patterns of the biases models-hindcast. The maps also show the mean biases averaged over the MS.   Table 2 reports the spatial correlations between the hindcast and the RCM data, computed through Pearson's coefficient.
The high values of the spatial correlation indicate that the spatial patterns are adequately captured by each model, for all the examined parameters and statistics. Altogether, these results show the ability of the models to reproduce reasonably well the wave dynamics on the historical period.

Analysis of trends on the future wave climate
First, the trends were computed on the ensemble time series of each wave parameter for each grid point. Trend analysis could be also performed directly on the series related to each single ensemble member, averaging the trend metrics at a second time. Such approach will be discussed in the second part of this section. Figure 2 summarizes the values of b related to H s and T m . First of all, it can be noticed how for all the examined statistics the vast majority of highlighted trends is negative, indicating an average reduction in the expected future wave parameters. According to the boxplots, positive values of b (i.e., upward trends) lie always above the upper bars, which denote the 75th percentile of the whole series. Moreover, the highest positive trends are generally less pronounced than the negative ones, especially as far as annual maxima data are concerned. Besides, the boxplots of Fig. 2 also reveal that the variability of b among the grid nodes increases with the order (mean, p90, and maxima) of the annual statistic considered. This can be clearly noticed by looking at the whiskers outside the upper and lower quartiles of the boxes, and leaves room for a further consideration: annual maxima data are more dispersed with respect to the annual lower percentiles, since the latter are influenced by mild sea states which, being more frequent, matter the most. This affects the estimations of b which show a higher variability for the annual maxima than for the other two statistics.
The aforementioned considerations are confirmed by the spatial distribution of b in the MS, as reported in Fig. 3. In this case, positive and negative values of b were linearly scaled in the [0,1] and [−1, 0] ranges, respectively, to better compare the trends location rather than the trends magnitude. To this end, we used the general formula: are located around the Balearic islands, in the southern part of the Ionian Sea and in the South-East Mediterranean basin around Cyprus. As for T m , the lowest values of b (i.e., the strongest downward trends) are found in the west coast of Greece and Crete. The weakest trends for annual mean and annual 90 th percentile series of both H s and T m are located in the Adriatic Sea, the North Tyrrhenian Sea, the Aegean Sea (but a small area characterized by positive trends, as previously highlighted) and east of the Tunisia coastline. When annual maxima are taken into account, the spatial distribution of expected trends becomes much more irregular (panels E and F for H max s and T max m , respectively). Positive trends are scattered over isolated spots in the Adriatic Sea, the Ionian Sea and off the south coast of Tunisia, while no broad areas uniformly characterized by significantly low values of b can be detected. Figure 3 allows assessing where trends characterized by different magnitude are most likely expected to take place, but it does not provide any information about their significance. To this end, the reliability of the b estimates was further evaluated looking at their 90% confidence interval, and coupling this information with the values of p MK computed on the respective time series, as explained in Sect. 2. Results related to the H s annual statistics are presented in Fig. 4, while results related to the T m annual statistics are presented in Fig. 5. The locations showing a change of sign between the upper and the lower confidence intervals of b are underlined in the leftmost side of the figures.
As regards H mean s and H max s , the areas characterized by non-significant trends are similarly located. In particular, attention is posed to the strait of Gibraltar and the Aegean Sea, being characterized by widespread marked areas (panels A and C of Fig. 4). A majority of the time series within such areas are in turn characterized by close-to-1 values of p MK , indicating the absence of a significant trend. A similar result characterizes the upward trends of H max s (panels E and F of Fig. 4). In fact, positive values of b in the Adriatic and in broad areas to the east of Tunisia and the Ionian Sea, are associated to change of sign between the respective b − and b + . A clear correlation exists also when negative trends are considered, that is, the values of b closest to 0 occur jointly with high values of p MK and no significant trends (e.g., in the southeast Mediterranean basin, the areas to the east of the Balearic islands, and in the North Tyrrhenian Sea).
As regards the time series of T m (Fig. 5), most of the positive trends in the Aegean Sea are still found to be not significant in case of annual mean (panels A and B) and annual 90th percentile (panels C and D). Given that the wave model was only forced by the RCM wind data, the trends previously highlighted are to be related to expected variations in the atmospheric circulation patterns. The MS is an enclosed basin at the mid-latitudes; thus, there is no environmental forcing affecting the wave climate as much as the wind (for instance, swells generated elsewhere or variation in the ice coverage). Overall, the results presented so far indicate a widespread decrease in the future wave heights and periods over the MS, and this seems to be consistent with expected decreases in the surface wind speed in this area, as shown for example by Mori et al. (2013) and Casas-Prat et al. (2018) (even though the latter considered a shorter period for the future projection, i.e., 2081-2100).
It can be clearly noticed that, when annual maxima are taken into account, the number of locations where trends are not significant increases dramatically (panels E and F in Figs. 4 and 5). This is due to the fact that the most negative trends are concentrated in isolated spots (as shown in panels E and F of Fig. 3); thus, it is difficult to detect homogeneous behaviors at basin level, compared with T mean m and T p90 m . Moreover, annual maxima refer to single instantaneous values per year, resulting in noisier datasets. Even though the methodology of Sen (1968) and Theil (1992) is sound with respect to possible outliers, this affects the computation of the confidence intervals, often leading to not significant trends.
In addition to the analysis so far described, as a compulsory step in the evaluation of the significance of the projected changes, we looked at the consistency of the trends related to the each single ensemble member. To this end, we computed the number of members resulting in concordant/discordant values of b for all the time series analyzed. Such number ranges from 7 − , meaning that all the models present negative trends, to 7 + (all the models presenting positive trends). Results are shown in Fig. 6.
A small stripe of fully consistent positive trends is found in the South Aegean Sea for H mean s and H p90 s , which exactly overlap with the locations characterized by significant upward trends according to panels A, B, C, and , a remarkable correlation can be noticed between the areas close to the Strait of Gibraltar, being characterized by significant positive trends according to both the approaches. Extending the analysis to T p90 m , a wide area east of Sardinia shows poor agreement among the ensemble members (panels C and D of Fig. 6), and this is exactly corresponding to the areas characterized by not significant trends underlined in panels A, B, C, and D of Fig. 5. Again, results related to the time series of annual maxima show a much more irregular spatial distribution (see panels E and F of Fig. 6).
As a further step, b was computed on the series of the 2006-2100 annual statistics averaged over the entire MS (i.e., across all the hindcast nodes taken into account), to derive an insight on the expected variation in the regional future wave climate. The ensemble values of regional b were carried out in 2 ways: (a) computing the slope of the ensemble values of the wave parameters. This estimate is hereinafter referred to as b r1 ; (b) computing the slope separately for each ensemble member, and then averaging on the ensemble members. This estimate is hereinafter referred to as b r2 . Figures 7, 8, and 9 show the time series of the regional H s statistics, while Figs. 10, 11, and 12 show the regional statistics relative to T m . In all the panels, the whole period covered by the models (1970-2100) is shown, along with the historical series computed from the hindcast developed and validated by Mentaschi et al. (2015). It can be seen from the orange lines in Figs. 7,8,9,10,11, and 12 that the hindcast series fall within the envelopes of the simulated ensembles for both H s and T m . Accordingly, the trends computed on the hindcast series are within the ranges provided by the estimates of b related to the ensemble members and summarized in Table 4 (see the Appendix), and serve to confirm the validity of the simulation models in producing reliable future wave climate projections. Results of b r1 and b r2 are reported in Table 3. From Figs. 7,8,9,10,11,and 12, it is clear at a glance that the regional time series of the investigated parameters are characterized by a downward trend. This finding can be also verified by the ITA, performed over the mean data of the regional averaged statistics (i.e., the black lines in Figs. 7,8,9,10,11,and 12)        show negative trends. As for the series of regional H max s , the ITA seems to be noisier, with a couple of data being above the no-trend line, though a clear negative trend can be pointed out also in this case.
Differences in the trends according to the ITA can be better appreciated by looking at the ecdf of δ, as shown in Fig. 15. Here, it can be noticed how the δ related to the annual maxima are generally farther from the zero line with respect to the other statistics, indicating a stronger negative trend. Similarly, the 90th percentile series are characterized by a trend more negative than the mean data. This particularly applies to the series of regional H s (panel A), while in case of T m the ITA trends are more similar among the investigated statistics (panel B). Such finding is further supported by the results of b r1 and b r2 , which show values increasing from the annual mean to the annual maxima.
These results show also that the differences between b r1 (trends computed on the parameters averaged across the models) and b r2 (trends computed as mean of the single model trends) are small. Indeed, the highest difference is observed for the time series of regional H max s , where using Finally, the mean direction θ mean m was considered. The trend analysis was carried out on the time series averaged over the whole MS, as shown in Fig. 16. Here, all the ensemble members agree on a slight clockwise shift.
A more detailed analysis was performed on six locations selected among different sub-basins of the MS, plotting the θ mean m series versus the respective years on a polar plot, along with the best linear fit. The investigated locations and the correspondent polar plots are shown in Figs. 17 and 18, respectively.
On average, an eastward trend seems to characterize the wave directions at locations A, D, and E, corresponding to the Gibraltar, South Mediterranean, and the Aegean basin, respectively. Conversely, no trends can be highlighted at location B (North Tyrrhenian Sea) and F (in front of the Nile Delta), while at location C scatters seem to be too dispersed to derive a reliable trend analysis, even though an anti-clockwise shift can be pointed out. This is most likely due to the local climatology of the Adriatic Sea, which is influenced by several local circulation patterns, resulting in no prevailing fetches for the point at hand.

Conclusion
In this work, we performed a trend analysis on time series of wave parameters projected up to the end of the twenty-first century in the MS. To this end, we used wave simulations driven by seven RCM over the 1970-2100 period under the RCP8.5 emission scenario. The future trends on time series of annual mean, annual 90th percentile, and annual maxima wave heights and periods were assessed through the Theil-Sen slope and the Mann-Kendall test, while the analysis of the wave directions relied on the use of polar plots. As for H s and T m , the results show that the trend metrics employed are generally consistent with each other, indicating the robustness of the projected changes. Such  Fig. 17 consideration is further confirmed by the similarities in the spatial patterns of the b estimates computed either on the parameters averaged across the models, and as the mean of the trends resulting from each single model.
Overall, the analysis revealed that the wave climate of the MS will be mainly characterized by downward trends, implying a progressive reduction in the magnitude of wave heights and periods. These results are in line with previous studies developed under the same emission scenario at larger scale, both for annual mean (Perez et al. 2015) and annual maxima (Wang et al. 2014). The trends computed on the time series of annual mean and annual 90th percentile are characterized by similar magnitudes and similar spatial distributions. On the other hand, the trends of annual maxima are more uncertain and irregularly distributed through the basin, despite the fact that the slope generally attains higher values. This suggests that the annual maxima should be used with caution, as they could result in dispersed time series, and lead to unreliable estimations of future trends.
Finally, as far as wave direction is concerned, a slight eastward trend is expected, but such behavior is not homogeneous across the different sub-basins.