Wave Analysis for Offshore Aquaculture Projects: A Case Study for the Eastern Mediterranean Sea

: The investigation of wave climate is of primary concern for the successful implementation of offshore aquaculture systems as waves can cause signiﬁcant loads on them. Up until now, site selection and design (or selection) of offshore cage system structures on extended sea areas do not seem to follow any speciﬁc guidelines. This paper presents a novel methodology for the identiﬁcation of favorable sites for offshore aquaculture development in an extended sea area based on two important technical factors: (i) the detailed characterization of the wave climate, and (ii) the water depth. Long-term statistics of the signiﬁcant wave height, peak wave period, and wave steepness are estimated on an annual and monthly temporal scale, along with variability measures. Extreme value analysis is applied to estimate the design values and associated return periods of the signiﬁcant wave height; structures should be designed based on this data, to avoid partial or total failure. The Eastern Mediterranean Sea is selected as a case study, and long-term time series of wave spectral parameters from the ERA5 dataset are utilized. Based on the obtained results, the most favorable areas for offshore aquaculture installations have been identiﬁed. cannot capture sufﬁciently due to the coarse spatial resolution. Comparing the obtained values of the evaluation metrics with the results from other studies (dealing with other pairs of wave datasets) focused on the Mediterranean Sea (e.g., [41]), it can be concluded that the present model performs better and provides more accurate values for H S and T P . From the obtained results, it is evident that model results tend to underestimate in situ measurements for both examined wave parameters.

1. Introduction 1.1. Aquaculture: Shifting from Coastal to Offshore Aquaculture production in the European Union (EU) represents only 1.0% and 1.5% of the global aquaculture production in terms of weight and value, respectively. Although EU fisheries production has shown a decreasing trend for the period 1990-2018, the aquaculture sector has gained ground in the seafood market, showing an increase of 24% from 1990 (up to 2018), but with a rather slow increase (6%) since 2007. The EU countries with the highest contribution in aquaculture production volume are Spain (27%), France (18%), Italy (12%), and Greece (11%), comprising 69% of the total EU production [1]. The Mediterranean Sea is an important contributor to this production, especially for seabream, seabass, and the Mediterranean mussel.
Fish are typically cultivated in cages located on land, in fresh water and sheltered coastal areas [2]. However, such protected environments tend to have limited carrying capacities, and there is a risk of harmful accumulation of farm waste products [3]. On the other hand, the future demand for fish is expected to increase. This is mainly attributed to the expected population increase, the gradual change in nutritional preferences towards the consumption of healthier food, and the reduction in the environmental footprint of the global food value chain, in particular the livestock industry [4].
Taking into consideration the above facts along with the scarce availability of coastal areas suitable for coastal aquaculture, it is logical to expand aquaculture production to offshore environments in order to reduce risks of disease for the stocks, provide sufficient dispersion of nutrients due to the presence of stronger currents [5], and mitigate the potential negative environmental consequences associated with coastal and land-based farms [6]. Moreover, operations on a larger scale are feasible, leading to an increase in annual production levels [7]. Let it be noted that exposed sites could be ideal for the cultivation of open-ocean species of high value, which can be supported by access to larger volumes of high-quality water [8]. In this framework, it is essential to make thorough research and adopt adaptive management for the sustainable expansion of offshore aquaculture [9,10].
On the other hand, the more intense environmental conditions, primarily attributed to the waves, ocean current and wind action, induce higher structural loads on fish cage structures and moorings, which can undergo failure if not properly designed. This, in turn, creates the need to design more robust fish cages, leading to more expensive solutions [11], while feeding, operation, and maintenance activities should be cautiously planned. As regards the influence of excessive wave action on the cultured fish, it is important to install fish cages at a water depth that will not be harmful to fish during storms and will not cause reduced growth and physiological problems under expected wave conditions. In addition, a priori knowledge of wild fish populations is important, because their interaction with cultured fish may cause several issues (e.g., spreading of parasites, genetic impacts) [6]. Other open issues include the limited management and monitoring capability of cages and the surrounding environment, water quality, and marine organisms that can be more cumbersome due to the distance from shore; however, some of the above can be mitigated with remote operations and monitoring, as well as automated feeding.
In order to identify whether a sea area represents a safe environment and is suitable for the development of offshore aquaculture installations, a detailed analysis of the wave climate is required as waves are considered one of the primary technical considerations [12,13]. Wave and current loads play a critical role for offshore fish cages, affecting both mooring loads and internal stresses. In a rapidly changing environment with short and steep wave patterns and more frequent manifestation of extreme wave events, it is imperative to examine seasonal variations in wave parameters and their corresponding extreme behavior, which is expected to greatly impact such structures, determine their operational success, and can lead to structure failures if not properly considered. The probabilistic design of offshore aquaculture systems can follow methods that are used in the design of coastal and offshore structures, which depends on the estimation of the return value of wave height [14][15][16][17].

Literature Review
Studies focused on the subject of site selection for offshore aquaculture are limited. In [18], a methodology was presented for the characterization of wave climate in the medium and long-term time scales to select the appropriate cage system and site of deployment around Tenerife Isl. (Canary Isl.). The analysis was based on 5-year time series of wave data from WAM model. The results were used to create a suitability map and identify the most favorable zones for the selected aquaculture cage systems. In [19], the most crucial criteria were identified for the site selection of open ocean aquaculture, extending from environmental parameters (e.g., hydrography, topography, meteorology) to socio-economic factors (e.g., manpower, infrastructure) and regulations. More recently, in [13], some comments were made on the challenges of offshore aquaculture along the Indian coastline, highlighting that the most critical physical parameters are wave height, current speed, bathymetry, and distance from shore. Other studies combined an analysis for site selection with a GIS software, e.g., [20][21][22][23][24][25]. Moreover, in [26], topographic, physical, chemical, biological, and oceanographic factors and various constraints were considered for the decision-making analysis along the Arabian Sea in Gujarat state (India), and in [27], a suitability index was proposed for finfish aquaculture within the Italian Economic Exclusive Zone. Finally, Zikra et al. [28] assessed environmental factors (including wave height) and water quality for fish cage site selection at Prigi Bay, Indonesia. It should be emphasized that after analyzing all the determinative factors in the context of a site identification study, a more detailed analysis must be carried out with in situ data of higher spatiotemporal resolution; see, e.g., [19]. However, this detailed site characterization phase is not discussed further in this paper.

Aim of Work
Although most of the studies already undertaken have considered wave action, only some descriptive statistics were provided at a local/regional level. In this study, wave analysis in the context of offshore aquaculture is taken a step further by recommending specific thresholds for various long-term statistics and design values. Another feature that distinguishes this work is the spatiotemporal coverage. The analysis is performed in a sub-basin level to detect promising areas for offshore aquaculture activities based on freely available long-term wave model data.
Specifically, the main aim of this study is to propose a methodology for the in-depth analysis and characterization of wave climate in an extended sea area in the context of the preliminary identification of favorable sites for the development of offshore aquaculture farms. The eastern Mediterranean Sea is selected as a case study, as the aquaculture sector in the corresponding EU countries (e.g., Italy, Greece, Cyprus) is well-organized and an opportunity exists to further increase the quantities of aquaculture production to meet the increasing demand for seafood [29]. As waves are highly variable in space and time, applications related to offshore mariculture usually require frequent wave measurements. The extended spatial context and the lack of in situ long-term high-resolution measurements implies the utilization of numerical model results. Therefore, spectral wave parameters of the ERA5 global reanalysis data have been utilized. This dataset was produced by the European Centre for Medium Range Weather Forecasting (ECMWF) and is composed by hourly wave parameters at 0.5 • horizontal resolution with a temporal coverage of 42 years (1979-2020), long enough to obtain as accurate as possible model-derived estimates. Apart from the standard parameters that characterize a wave climate in a region (i.e., significant wave height H S , peak or mean wave period T P , T m ), spectral wave steepness s m is also considered in this study. The behavior of offshore structures is highly dependent on the wave steepness; as spectral wave steepness increases, the probability of experiencing instabilities in the sea surface gets higher [30]. Moreover, rogue (or freak) waves can be triggered by a sea state of high steepness, but additional factors should also be present for a higher probability of occurrence of such waves, like a narrow-banded wave spectrum in frequency and direction [31]. Consequently, the design and choice of offshore fish cages are directly related to the wave steepness.
The remaining part of this paper is structured as follows: Section 2 describes in detail the model data used in this study and their validation against buoy measurements. Section 3 presents the theoretical background of the wave parameters used in this study, along with the basic concepts of Block Maxima in extreme value theory. Moreover, the proposed methodology for the identification of favorable sites for the development of offshore aquaculture farms is described in the last subsection. Section 4 discusses the obtained results from the wave climate analysis and the proposed methodology. Some additional aspects are discussed in Section 5; conclusions and future considerations are presented in Section 6.

Description
The latest global wave reanalysis dataset, produced by ECMWF, ERA5, was obtained from the Copernicus Climate Data Store portal (https://cds.climate.copernicus. eu/#!/home, accessed on 20 September 2021) in the form of a regular grid continuously over time and space, from January 1979 until December 2020 [32]. Hourly data are provided for this dataset at a 0.5 • spatial resolution for the two horizontal dimensions. The dataset is based on the Integrated Forecasting System (IFS) Cy41r2 and uses a 4D variational data assimilation scheme of satellite and in situ observations, which are substantially more than its predecessor, ERA-Interim; for more details, see also [33]. The considered spectral wave parameters are the significant wave height H S , the peak wave period T P , and the spectral wave steepness s m . Although the spatial resolution of this dataset is coarse, it is suitable for the purposes of this study, which is focused on a rather extended offshore sea area. Our approach is not expected to be useful near the coast (i.e., lower than 10 km) with this dataset, but could be applied to higher-resolution data for that purpose.
The validity of ERA5 wave data has been studied in many works for different parts of the world. For instance, Mahmoodi et al. [34] evaluated ERA5 wave data with two buoy stations in the Persian Gulf by using four goodness-of-fit metrics and concluded that the dataset provided precise results. A good agreement between the measured buoy data and the ERA5 results for significant and maximum wave heights (bias around 0.29 m and 0.18 m for coastal and deep waters, respectively) was also reported by Muhammed Naseef and Sanil Kumar [35] for the Bay of Bengal. In [36], significant wave heights and mean wave periods from the ERA5 dataset were evaluated based on four statistical parameters at one buoy located off the Oman coast. From that study it was revealed that the ERA5 swell wave heights were overestimated during monsoon and post-monsoon seasons, while wind wave heights were not accurate for the monsoon period (RMSE = 0.32 m). On the other hand, Hisaki [37] evaluated ERA5 wave height data with wave measurements from buoys at 18 coastal locations around Japan and concluded that ERA5 wave height tends to be underestimated for increasing wave heights.

Evaluation
The evaluation and validation of the ERA5 dataset at the examined sea area is based on in situ measurements from 10 oceanographic buoys obtained from the POSEIDON marine monitoring network of Greece, the Italian national data buoy network (RON -Rete Ondametrica Nazionale) and the closest grid points of the ERA5 results. The time series of the wave measurements (for H S and T P ) were obtained from the Copernicus Marine Service portal (https://resources.marine.copernicus.eu/products, accessed on 27 September 2021) with the product identifier INSITU_MED_NRT_OBSERVATIONS_013_035. The location, depth, distance from the closest ERA5 grid point, and period of measurements are shown in Table 1; see also Figure 1.    Buoy measurements characterized as 'good' in the quality control were kept, while stuck values were disregarded from the analysis. Moreover, outliers were identified according to the studentized residuals, which take into account deviations in both variables [38], and ignored after thorough examination of the histogram of the studentized residuals and scatter plots between measured and modeled variables. The maximum percentage of outliers that was detected among the examined data samples was 2.5%. Despite the above filtering and time collocation of the two data sources, the remaining number of measurements was sufficient to carry out the statistical analysis. In Tables 2 and 3, the main descriptive statistics (sample size N, mean value m, standard deviation s, minimum min, maximum max, coefficient of variation CV, skewness Sk, and kurtosis Ku) of H S and T P are presented, respectively, for both data sources at the examined locations. Considering all the locations, the mean values of H S seem to be underestimated by the ERA5 wave model, apart from SAN, while the overall max values of H S are higher for the buoy measurements, apart from ANC and ATH. CV values are higher for the model except for CRO, MAZ, MYK, and SAN; on the other hand, the values of s are higher for in situ measurements apart from SAN and ZAK. The positive values of Sk higher than 1 denote a highly skewed (to the right) distribution for H S and the relatively high values of Ku, especially for CRO, indicate heavy tails on the right side.
As regards T P , the mean values are lower for the model data, apart from MYK and SAN; max values, as well, are lower for all locations. CV values are rather close for the two data sources (apart from MYK and SAN). Sk values are all positive (except for model data at SAN) and close to 0.5, indicating a fairly symmetrical distribution, which is also supported by the values of Ku around 3.
The goodness-of-fit metrics used for the evaluation are the following: correlation coefficient (R), mean error (BIAS), root mean square error (RMSE), mean absolute error (MAE), scatter index (SI), and symmetric mean absolute percentage error (SMAPE). For their definition, see [39,40]. The model performance is characterized as good if all the values of the above metrics are as low as possible, apart from correlation coefficient that should be close to unity. In Figures 2 and 3, the scatter density plots for H S and T P , respectively, are presented for all examined locations, along with the values of the evaluation metrics. Regarding H S , R values are above 0.91 for all locations, apart from SAN, indicating a strong positive linear relationship between the model and buoy data. BIAS is consistently positive (apart from SAN), verifying once more that this variable is underestimated by the ERA5 wave model. MAZ and ZAK for the Italian and Greek buoy network, respectively, present the lowest values for RMSE, MAE, SI, and SMAPE. Regarding T P , from the values of BIAS it is clear that T P is underestimated at all locations apart from MYK and SAN. ATH presents the lowest values for RMSE, MAE, SI, and SMAPE and the highest value for R. A possible explanation for the poor performance of the model at SAN is the rather large distance between the buoy and the closest grid point (26.4 km) in combination with the complex topography of the wider area (Cycladic Islands Complex), which the model cannot capture sufficiently due to the coarse spatial resolution. Comparing the obtained values of the evaluation metrics with the results from other studies (dealing with other pairs of wave datasets) focused on the Mediterranean Sea (e.g., [41]), it can be concluded that the present model performs better and provides more accurate values for H S and T P . From the obtained results, it is evident that model results tend to underestimate in situ measurements for both examined wave parameters.
network, respectively, present the lowest values for RMSE, MAE, SI, and SMAPE. Regarding , from the values of BIAS it is clear that is underestimated at all locations apart from MYK and SAN. ATH presents the lowest values for RMSE, MAE, SI, and SMAPE and the highest value for R. A possible explanation for the poor performance of the model at SAN is the rather large distance between the buoy and the closest grid point (26.4 km) in combination with the complex topography of the wider area (Cycladic Islands Complex), which the model cannot capture sufficiently due to the coarse spatial resolution. Comparing the obtained values of the evaluation metrics with the results from other studies (dealing with other pairs of wave datasets) focused on the Mediterranean Sea (e.g., [41]), it can be concluded that the present model performs better and provides more accurate values for and . From the obtained results, it is evident that model results tend to underestimate in situ measurements for both examined wave parameters.

Brief Theoretical Background for Wave Parameters
In deep water conditions, the spectral wave steepness in a sea state is defined by the ratio ⁄ , where is the wavelength; this, in turn, is related to wave period based on the dispersion relation. In irregular sea states and in terms of and (mean zerocrossing wave period), is given by the following relation: where = 4 and = 2 ⁄ , with = d , = 0, 1, 2, …, the spectral moments, is the single-sided frequency wave spectrum, is the wave frequency, and = 9.81 m/s 2 is the gravitational acceleration. Subsequently, can be calculated for the examined basin in terms of the hourly results of and . Other studies where spectral wave steepness was included for ocean engineering applications include [42,43]. Peak wave period corresponds to the spectral peak frequency, i.e., = 1 ⁄ .

Brief Theoretical Background for Wave Parameters
In deep water conditions, the spectral wave steepness s m in a sea state is defined by the ratio H S /λ, where λ is the wavelength; this, in turn, is related to wave period based on the dispersion relation. In irregular sea states and in terms of H S and T Z (mean zero-crossing wave period), s m is given by the following relation: . . , the spectral moments, S( f ) is the single-sided frequency wave spectrum, f is the wave frequency, and g = 9.81 m/s 2 is the gravitational acceleration. Subsequently, s m can be calculated for the examined basin in terms of the hourly results of H S and T Z . Other studies where spectral wave steepness was included for ocean engineering applications include [42,43]. Peak wave period T P corresponds to the spectral peak frequency, i.e., T P = 1/ f P .

Theoretical Background for Block Maxima Method in Extreme Value Theory
In the context of extreme value analysis, there are two commonly used approaches for estimating return period values for extreme environmental conditions and design purposes: the block maxima (BM) and the peaks-over-threshold (POT) methods. In the former approach, equal-sized non-overlapping bins (blocks) are generated to extract maximum observations from the data sample, while in the latter one, observations above a certain threshold, appropriately selected, are extracted. Both approaches are based on an independent and identically distributed (iid) assumption as regards the underlying sample (typically a stationary time series data set). The BM method is considered in this work as (i) the available dataset extends to several decades, ensuring a large data sample for curve fitting in the estimation of the 50-year return period, and (ii) there is not yet an established automated and robust methodology for the threshold selection of the POT method (see, e.g., the review in [44]).
Let M n be the maximum of a sequence of iid random variables X 1 , X 2 , . . . , X n , i.e., M n = max{X 1 , X 2 , . . . , X n }. From the extremal types theorem, (see e.g., [45]) if appropriate sequences of normalizing constants a n > 0, b n exist as n → ∞, then M * n = (M n − µ n )/σ n , follows the limiting distribution G(x): where G(x) is a non-degenerate function and belongs to one of three possible families of extreme value distributions; see also [46,47]. The three distributions can be unified in a single family, namely the generalized extreme value (GEV) distribution [48], expressed by the following cumulative distribution function: defined on the set {x : 1 + ξ(x − µ)/σ > 0} while the corresponding probability density function is where µ ∈ R, σ > 0, ξ ∈ R are the location, scale, and shape parameters, respectively. The shape parameter ξ governs the tail behavior of the distribution at its upper bound. For different values of ξ, this distribution is modified corresponding to different distribution types. Specifically, (i) for ξ → 0 , the Gumbel distribution is represented (type-I); (ii) for ξ > 0, the Fréchet distribution is represented (type-II); and (iii) for ξ < 0, the reversed Weibull distribution is represented (type-III). The estimation of the three unknown parameters of the GEV distribution can be relied on the maximum likelihood method and method of moments, among other things; see also [49,50], where a thorough review on the assessment of various methods on the estimation of the GEV parameters and the effects of the available sample size are presented. In this work, the L-moments method [51] was implemented for the estimation of GEV parameters, as there were some cases where the maximum likelihood estimators were not obtainable. For this purpose, the 'extRemes' package, available at the Comprehensive R Archive Network (https://cran.r-project.org/web/packages/extRemes/ accessed on 30 September 2021), was used.
Having a long sequence of independent observations X 1 , X 2 , . . ., and dividing it into equal blocks of length n, a series of block maxima M n,1 , M n,2 , . . . , M n,m , can be obtained and the GEV distribution can be fitted in this sample. The total number of observations is n × m.
Another concept in extreme value analysis is the estimation of the likelihood of rare events, e.g., the probability of exceeding, once in the mean, a threshold larger than the highest observation of M n . This probability is expressed through the return level x p associated with a specific return period T R , also known as the recurrence interval. The return period is the average time of occurrence of events of this level or greater, defined by The return level x p is estimated by the following equation: which actually corresponds to the (1 − 1/T R )− quantile of F. In terms of the GEV parameters (i.e., F is approximated by a GEV distribution G), x p is expressed as follows:

Methodology for the Identification of Favourite Sites for Offshore Aquaculture
According to the definition of offshore aquaculture proposed by [52], aquaculture is considered 'offshore' if it takes place in the open sea and has significant exposure to wave action, among other things, and 'where there is a requirement for equipment and servicing vessels to survive and operate in severe sea conditions from time to time'. Hence, accurate knowledge of wave conditions is necessary not only for the survival of floating cages and the fish, but for the humans and devices/vessels that support and supplement their operation and maintenance. To this end, the criteria that should be met as regards the wave conditions for the development of offshore aquaculture should consider all the above aspects.
Wave conditions and bathymetry belong to the so-called technical (or engineering) constraints. Other parameters related to this category are current speed and direction, as well as the type of seabed substrate, which are not utilized here as it is beyond the scope of this work. To the authors' knowledge, the only internationally accepted standard reference for the classification of the prevailing environmental conditions for marine fish farms was provided by [53], which has recently replaced [54]. However, access to the latest Norwegian regulations was not possible at the time of writing, and potentially useful information cannot be included herein. In [55], where operational limits were defined for aquaculture operations from a risk perspective, nine categories for significant wave height were suggested based on the World Meteorological Organization's codes for sea state. Since there is very limited data and information available in the literature on the maximum allowable values of wave parameters such aquaculture systems can withstand, the thresholds of the considered wave characteristics were set primarily considering the resistance of offshore aquaculture farms in a high-energy environment and the physics of ocean waves. Moreover, the limits for H S , provided in Table 4, are fairly well-aligned with the ones provided in [56], where H S ∈ [0.0, 0.5] corresponds to low exposure, H S ∈ [0.5, 1.0] corresponds to moderate exposure, and H S ∈ [2.0, 3.0] corresponds to high exposure; see also [18].
The considered wave parameters and water depth were classified in two categories; (i) optimal, where the proposed conditions are the most suitable ones for the cultivation of a variety of aquaculture species, and (ii) suboptimal, where the proposed conditions are suitable but could result in reduced productivity. The proposed thresholds for each parameter are summarized in Table 4. Note that the values for H S , T P , and s m refer to the annual time scale values for the identification of potential sites for further in-depth analysis. Supplementary and more accurate information on the wave conditions is given by variability measures and relevant analysis at finer time scales (e.g., monthly values). One of the most commonly used fish cages is the circular 'PolarCirkel' type, which can withstand significant wave heights between 3.5 m and 4.5 m, depending on the construction details and overall size. On the other hand, extreme wave heights can be a major risk factor and influence the development of offshore aquaculture. Given the relatively mild wave climate of the Mediterranean Sea, the threshold for H S for the optimal case was considered below 1 m, while for the suboptimal case, it was selected between 1 m and 3 m. Regarding the 50-year return period H S , for the optimal case, it was set up below 4.5 m to include 'PolarCirkel' types of fish cages and for the suboptimal case, a range between 4.5 m and 7 m was considered. The threshold value of 5 s was selected in the optimal case for spectral peak wave period, as values lower than this threshold are not likely to generate significant hydrodynamic loads on the aquaculture system as they are associated with low H S values (usually below 1 m in the examined area). From a numerical simulation study on the dynamic response of offshore fish cage systems, it was found that the mooring line tension was small for values of s m below 1/30 under regular wave loading [57]. As the wave breaking limit (s m > 0.14) is not likely to happen in deep water, steep waves with relatively high significant wave heights and short associated wave periods might occur. In this context, the optimal case for s m included values below 0.035 and the suboptimal case was taken in the range (0.035, 0.08).
In this study, water depth is considered between 50 m and 300 m; the lower threshold has been set by [58] based on the opinions of experts from a technical workshop and the upper one has been set by the existing technology of offshore fish cages (e.g., Ocean Farming AS). The recommendations in [13] were also considered for the upper limit of the optimal range for the depth profile. The bathymetric grid data for the examined basin were obtained by the EMODnet-Bathymetry portal and the 2016 EMODnet Digital Terrain Model (DTM) product with a resolution (for x-, y-direction) equal to 0.0021 • [59]. The reason for selecting this version, and not the most recent one, is to have a similar resolution to the spectral wave data.
The procedure of site identification for offshore aquaculture based on the time series of spectral wave parameters and gridded bathymetric data can be summarized in the following steps; see also Figure 4: Step 1: Estimation of the long-term average characteristics of the sea state. It is recommended to analyze H S , T P , and s m time series. The time scale can include annual, seasonal, and monthly values.
Step 2: Estimation of extreme wave parameters like the 50-year return period H S , i.e., H S value that, on average, is exceeded once every 50 years. Such values are crucial for the design of offshore aquaculture systems.
Step 3: Estimation of long-term variability measures of the above spectral wave parameters to quantify variations over time. Such measures typically include standard deviation, mean annual variability, and inter-annual variability. The definition of the two latter measures can be found in [60]. Areas with low temporal variability are preferred, since they are favorable regarding logistics, maintenance activities, scheduling of in situ visits, etc.
Step 4: Application of the thresholds presented in Table 4 in the eastern Mediterranean Sea to identify optimal and suboptimal locations for potential offshore aquaculture development.
Step 5: Investigation of the behavior, variability and statistics of finer time scales at the identified optimal and suboptimal locations.
Step 6: Final selection of optimal and suboptimal locations.
Climate 2022, 10, x FOR PEER REVIEW 13 of 26 Step 4: Application of the thresholds presented in Table 4 in the eastern Mediterranean Sea to identify optimal and suboptimal locations for potential offshore aquaculture development.
Step 5: Investigation of the behavior, variability and statistics of finer time scales at the identified optimal and suboptimal locations.
Step 6: Final selection of optimal and suboptimal locations.

Annual and Hourly Analysis
In Figure 5, the spatial distribution of various descriptive statistics for is presented. High mean annual values of (above 1.1 m) seem to coincide with high standard deviation values for the S Aegean Sea and the western part of Sicily. The Adriatic Sea and the E Levantine Sea are characterized by low mean annual and 90th percentile values of (below 0.6 m and 1.4 m, respectively). The E Levantine Sea also presents relatively low variability values. High mean annual variability (MAV) values of are depicted at the NW Aegean Sea and along the Italian coasts, while the relevant values for inter-annual variability (IAV) are located in the E Ionian and Tyrrhenian seas.  The corresponding descriptive statistics for are shown in Figure 6. Specifically, the Adriatic, Aegean, and NE Levantine seas present low mean annual and 90th percentile values compared to the southern part of the examined sea area. Low mean annual and inter-annual values are found in the S Aegean Sea and W Levantine Sea. The corresponding descriptive statistics for T P are shown in Figure 6. Specifically, the Adriatic, Aegean, and NE Levantine seas present low mean annual and 90th percentile values compared to the southern part of the examined sea area. Low mean annual and inter-annual values are found in the S Aegean Sea and W Levantine Sea. In Figure 7, the descriptive statistics for wave steepness are presented. Sea areas that combine low values of mean annual and low mean annual and inter-annual variability are the N Adriatic Sea, off the southern part of Crete, and Cyprus. Given the spatial distribution of the mean annual values and 90th percentile of , it is expected that sea states with high steepness (above 0.045 and 0.065, respectively) are more likely to occur in In Figure 7, the descriptive statistics for wave steepness s m are presented. Sea areas that combine low values of mean annual s m and low mean annual and inter-annual variability are the N Adriatic Sea, off the southern part of Crete, and Cyprus. Given the spatial distribution of the mean annual values and 90th percentile of s m , it is expected that sea states with high steepness (above 0.045 and 0.065, respectively) are more likely to occur in the Aegean Sea.

Monthly Analysis
In Figure 8, the spatial distribution of the mean monthly values of is presented. The highest mean monthly values of are observed during winter (up to 1.8 m for January and 1.7 m for December and February) off the southern part of the Ionian Sea. During this period, the entire examined area, apart from the Adriatic and NE Levantine seas, seems to be characterized by high-energy sea waves. During July, August, and September, the effect of the Etesian wind is evident at the SE part of the Aegean Sea (located at the strait between Crete and Rhodes), accounting for relatively high values (around 1.

Monthly Analysis
In Figure 8, the spatial distribution of the mean monthly values of H S is presented. The highest mean monthly values of H S are observed during winter (up to 1.8 m for January and 1.7 m for December and February) off the southern part of the Ionian Sea. During this period, the entire examined area, apart from the Adriatic and NE Levantine seas, seems to be characterized by high-energy sea waves. During July, August, and September, the effect of the Etesian wind is evident at the SE part of the Aegean Sea (located at the strait between Crete and Rhodes), accounting for relatively high H S values (around 1.3 m) compared to the rest of the eastern Mediterranean Sea. The western part of the examined area is characterized by low H S values (below 1 m) for five consecutive months (between May and September). The spatial distribution of mean H S for March and October is rather homogeneous, with higher values (up to 1.4 m) observed during March.
Climate 2022, 10, x FOR PEER REVIEW 1 Figure 8. Spatial distribution of mean for each month.
In Figure 9, the spatial distribution of the mean monthly values of is pres The highest mean monthly values of are observed during winter (around 7 s fo cember, January, and February) located off the eastern coasts of Libya. Novembe March follow with mean values around 6.5 s. April, May, and October seem to a homogeneous spatial distribution of mean , excluding the Adriatic and norther gean seas. The spatial distribution of mean is observed to have a similar pattern ing the summer months (June, July, August) and September: a gradual increase of from the NW to the SE direction. In Figure 9, the spatial distribution of the mean monthly values of T P is presented. The highest mean monthly values of T P are observed during winter (around 7 s for December, January, and February) located off the eastern coasts of Libya. November and March follow with mean T P values around 6.5 s. April, May, and October seem to have a homogeneous spatial distribution of mean T P , excluding the Adriatic and northern Aegean seas. The spatial distribution of mean T P is observed to have a similar pattern during the summer months (June, July, August) and September: a gradual increase of mean T P from the NW to the SE direction.
Climate 2022, 10, x FOR PEER REVIEW 18 Figure 9. Spatial distribution of mean for each month.
In Figure 10, the spatial distribution of the mean monthly values of is prese Generally, the Aegean Sea is characterized by relatively high values compared to the sea areas during all months. The highest values (0.056) are observed during July and gust, followed by winter (0.054). The Ionian, Adriatic, Tyrrhenian, and eastern Leva seas are areas with low mean values of for all months, especially during summe Figure 9. Spatial distribution of mean T P for each month.
In Figure 10, the spatial distribution of the mean monthly values of s m is presented. Generally, the Aegean Sea is characterized by relatively high values compared to the other sea areas during all months. The highest values (0.056) are observed during July and August, followed by winter (0.054). The Ionian, Adriatic, Tyrrhenian, and eastern Levantine seas are areas with low mean values of s m for all months, especially during summer.

Extreme Value Analysis of Significant Wave Heights
Based on Figure 11a, the examined basin is characterized overall by values of the location parameter for GEV distribution above 4, apart from the Adriatic Sea and some coastal areas around Cyprus, Greece, Italy, Tunisia, and Turkey, where extreme H S demonstrates milder behavior; thus, at a first glance, the above-mentioned areas seem to be more suitable for offshore aquaculture development. From Figure 11b, the spatial distribution of the scale parameter of the GEV distribution denotes higher variability in the yearly extreme H S values mainly in the northern Aegean and eastern Tyrrhenian seas, and off the western part of Libya. Due to this behavior, these areas could be excluded from an in-depth analysis. The most important parameter of GEV distribution in terms of very rare events is the shape parameter ξ, since it determines the right-tail behavior of the distribution; therefore, the sign of ξ is of importance when extrapolating to long return periods. Based on Figure 11c, the spatial distribution of ξ is not homogeneous. Specifically, most of the grid points located in the examined sea area have a negative value for the shape parameter (up to −0.6), corresponding to a reversed Weibull distribution. The Aegean Sea, excluding its middle part, and the SW Levantine Sea are characterized by the highest absolute values of the shape parameter for GEV distribution, indicating that the quantiles of the H S distribution with a negative shape parameter tend to a finite upper bound. Figure 11d, showing the 50-year design values of H S , presents a similar pattern to the location parameter for GEV distribution. The highest value (8.6 m) is located in the Tyrrhenian Sea, while other areas with values over 8 m are below the 36th parallel, between the 13th and 22nd meridians. Apparently, the abovementioned sea areas are not favorable for the deployment of offshore aquaculture systems. Mild exposure to extreme wave heights is observed for the Adriatic Sea, the northern Aegean Sea, some nearshore areas in the NE Levantine Sea and off the coast of Tunisia.
Climate 2022, 10, x FOR PEER REVIEW 1 Figure 10. Spatial distribution of mean for each month.

Extreme Value Analysis of Significant Wave Heights
Based on Figure 11a, the examined basin is characterized overall by values location parameter for GEV distribution above 4, apart from the Adriatic Sea and coastal areas around Cyprus, Greece, Italy, Tunisia, and Turkey, where extrem demonstrates milder behavior; thus, at a first glance, the above-mentioned areas se be more suitable for offshore aquaculture development. From Figure 11b, the spati tribution of the scale parameter of the GEV distribution denotes higher variability  favorable for the deployment of offshore aquaculture systems. Mild exposure to extreme wave heights is observed for the Adriatic Sea, the northern Aegean Sea, some nearshore areas in the NE Levantine Sea and off the coast of Tunisia.
(a) (b) (c) (d) Figure 11. Spatial distribution of (a) location parameter, (b) scale parameter, (c) shape parameter of the GEV distribution, and (d) 50-year design value for .

Identification of Sites Based on Wave Climate and Bathymetry
After applying the thresholds presented in Table 4, two sites were identified as optimal for further investigation, which are in the NE Adriatic and N Ionian seas, while five locations were found to meet the specified requirements for the suboptimal case, all located in the southern part of the eastern Mediterranean Sea; see also Figure 12. Let it be remembered that for the selection of the sites, the bottom depths have also been considered. Taking into consideration the results of the variability measures and the mean monthly values for at the identified sites, shown in Figure 13, the following remarks can be made: • For the optimal sites, O1 site is characterized by a very mild wave climate for all months due to its sheltered location, with values below 0.41 m. For both sites, MAV

Identification of Sites Based on Wave Climate and Bathymetry
After applying the thresholds presented in Table 4, two sites were identified as optimal for further investigation, which are in the NE Adriatic and N Ionian seas, while five locations were found to meet the specified requirements for the suboptimal case, all located in the southern part of the eastern Mediterranean Sea; see also Figure 12. Let it be remembered that for the selection of the sites, the bottom depths have also been considered. Taking into consideration the results of the variability measures and the mean monthly values for H S at the identified sites, shown in Figure 13, the following remarks can be made: Similar results are drawn for the other two spectral wave parameters, and (figures are not presented), in terms of the optimal sites.  optimal locations, with S2 and S3 having the lowest IAV and MAV values, respectively.
Similar results are drawn for the other two spectral wave parameters, and (figures are not presented), in terms of the optimal sites.  Similar results are drawn for the other two spectral wave parameters, T P and s m (figures are not presented), in terms of the optimal sites.

Discussion
According to this study, seven locations were characterized as preferable to be further investigated for the potential development of offshore aquaculture based on wave variables and bathymetry. The small number of locations identified by the methodology and their sparse spatial distribution is attributed to the coarse spatial resolution of the wave model grid and the steep continental shelf of the Mediterranean Sea. There are many sea areas where the model does not cover the isobath zone extending from 50 m to 300 m that was selected for the purposes of this study. Therefore, the water depth criterion eliminates many sea areas that have suitable wave conditions. Another important and well-known issue is the uncertainty associated both with the wave data source (i.e., modelling uncertainties) and the extreme value analysis (i.e., statistical uncertainties). As was mentioned in Section 2.2, wave model results seem to underestimate buoy measurements for both H S and T P , and tend to underpredict large waves and extreme n− year significant wave heights. This behavior could be attributed to the limitations of most wind reanalysis datasets that are not able to resolve high-energy wind fields of a finer spatial scale, e.g., [61]. On the other hand, the statistical uncertainties include the application of the fitting method (both model and block size) used in this study, i.e., the BM method with block sizes corresponding to one year, which drastically reduces the sample size, e.g. [62]. For the sake of completeness, the estimation of the 95% confidence intervals (CIs) of the 50-year return values H S was performed for the optimal and suboptimal locations. The parametric bootstrap was used by generating 1,000 bootstrap samples and the results are presented in Table 5. The upper 95% CI limit falls under the threshold specified for the 50-year return period H S at the optimal locations. For the suboptimal sites, the upper 95% CI limit exceeds the defined relevant threshold by 0.65 m at the highest. Considering the above uncertainties and given that the deterministic design of offshore aquaculture systems follows methods similar to coastal and offshore structures, a safety margin can be incorporated in this case, which is represented by multiplying the 50-year design values of H S with an appropriate safety factor (usually between 1.1 and 1.3) [17].

Conclusions
The present work proposed a methodology for the identification of potential sites to develop offshore aquaculture systems in the eastern Mediterranean Sea. The identification was based on two key technical factors: (i) the characterization of the wave climate, and (ii) the water depth. To this end, a long-term wave dataset was used, the ERA5 wave model, with hourly estimates of significant wave height and peak wave period, which are fundamental for the structural behavior of offshore fish cages. Wave steepness was also included in the study as it influences the design or selection of offshore fish cages. The validity of the model results was ensured by thoroughly assessing the performance of the wave model with in situ measurements obtained from 10 offshore and coastal buoys in the Italian and Greek seas, in terms of significant wave height and peak wave period. The wave analysis was carried out in annual and monthly time scales to reveal different patterns in terms of both long-term statistics and variability measures under normal operating conditions. Moreover, the 50-year significant wave height was used to characterize extreme wave conditions under the assumption of stationary wave climate trends. Threshold values were recommended for the above wave statistics by taking into consideration current offshore aquaculture technologies and the physical characteristics of waves. Finally, optimal and suboptimal locations were identified in the Adriatic and Ionian seas, and along the Mediterranean African coastline, respectively.
Based on the results of this study, the following findings are summarized: • The thorough evaluation of the ERA5 wave model data showed good agreement with buoy measurements in the basin of interest. • Wave analysis in extended sea areas for site selection purposes should be based on long-term data and include average statistics, variability measures and extreme values.
• Apart from common sea state variables, wave steepness is also recommended to be embedded in relevant studies as an additional variable affecting the design of offshore fish cages.

•
The proposed methodology is a step forward in supporting the sustainable development of offshore aquaculture in the Mediterranean Sea and elsewhere. It can also be considered as a precursor of a more holistic approach, incorporating more parameters that influence the site selection procedure.
It is evident that the proposed methodology can be easily applied to any extended sea area (e.g., the Black Sea and the West Mediterranean Sea) on a national/trans-national level in the context of a pre-evaluation site selection process for offshore aquaculture development. Moreover, it could be applied on a regional or local level after appropriate modifications. For instance, in the latter case, high-resolution long-term data sets are required as they provide more accurate estimates near the coast. In case of more limited sea areas, calibration of the model results can also be performed to improve the accuracy of the statistical measures, given that the buoys are in approximately the same types of environments under investigation for offshore aquaculture. It is highlighted that the recommended thresholds for the spectral wave parameters should be used with caution for high-energy sea areas.
In the future, it would be valuable to include additional criteria, such as currents, distance from shore, proximity to ports, ship routes and other marine uses, and implement a multi-criterion decision-making analysis to identify priority sites for the development of offshore aquaculture projects. Such results could highlight possible obstacles and assist policy makers in promoting offshore aquaculture through a holistic approach. Moreover, a downscaled approach with additional criteria will allow the in-depth examination of site-specific factors that will determine additional aspects of offshore aquaculture activities, such as carrying capacity and cost-benefit findings.