A Raster-Based Multi-Objective Spatial Optimization Framework for Offshore Wind Farm Site-Prospecting ^†

Abstract

Siting an offshore wind project is considered a complex planning problem with multiple interrelated objectives and constraints. Hence, compactness and contiguity are indispensable properties in spatial modeling for Renewable Energy Sources (RES) planning processes. The proposed methodology demonstrates the development of a raster-based spatial optimization model for future Offshore Wind Farm (OWF) multi-objective site-prospecting in terms of the simulated Annual Energy Production (AEP), Wind Power Variability (WPV) and the Depth Profile (DP) towards an integer mathematical programming approach. Geographic Information Systems (GIS), statistical modeling, and spatial optimization techniques are fused as a unified framework that allows exploring rigorously and systematically multiple alternatives for OWF planning. The stochastic generation scheme uses a Generalized Hurst-Kolmogorov (GHK) process embedded in a Symmetric-Moving-Average (SMA) model, which is used for the simulation of a wind process, as extracted from the UERRA (MESCAN-SURFEX) reanalysis data. The generated AEP and WPV, along with the bathymetry raster surfaces, are then transferred into the multi-objective spatial optimization algorithm via the Gurobi optimizer. Using a weighted spatial optimization approach, considering and guaranteeing compactness and continuity of the optimal solutions, the final optimal areas (clusters) are extracted for the North and Central Aegean Sea. The optimal OWF clusters, show increased AEP and minimum WPV, particularly across offshore areas from the North-East Aegean (around Lemnos Island) to the Central Aegean Sea (Cyclades Islands). All areas have a Hurst parameter in the range of 0.55–0.63, indicating greater long-term positive autocorrelation in specific areas of the North Aegean Sea.

1. Introduction

Offshore wind energy (OWE) was first developed in the early 1990s in the North European countries and, during the last three decades, has expanded with tremendous efforts. Development has been further accelerated by 2000 and on, when larger wind turbines became available, water depth and distance from the shore increased, and experience was gained with the first operational OWFs [1]. In the early stages of an offshore wind farm (OWF) deployment, the industry and policymakers considered it a straightforward alternative to the challenges of increasing onshore wind capacities [2]. However, the outcomes from offshore wind projects were often less favorable than anticipated. The extended planning periods and the rising costs of manufacturing, installation, and operations attracted significant criticism [3]. Despite these challenges, research on offshore wind energy (OWE) has deepened our understanding of the factors influencing offshore wind investments. Moreover, the complex and interdisciplinary nature of OWF site-prospecting campaigns has necessitated a thorough understanding and careful consideration of the criteria for selecting suitable marine areas.

Focusing on the previously mentioned characteristics, OWE is characterized by a spatio-temporal nature, where changes are obtained at many different scales, including energy production and available marine space. The spatial limitations are absolute regarding areas excluded from development because of legal, technical, or physical reasons (wind, depth, or other marine activities) [4]. This indicates that a hierarchical selection must be established for a successive planning scheme, focusing on the sites’ optimality for future OWE deployment. With technological improvements and the increased energy demand, areas considered less attractive during the past years are subject to future development. Accordingly, the temporal perspectives of wind speed and wind power potential and variability, combined with the economics of OWE, mean that appropriate methods are required for monitoring and quantifying many critical aspects, which operate continuously in time and space [3,5].

Thus, the offshore wind sector is considered a complex decision environment, where decisions must be taken on both the strategic and operational levels at various stages of the life cycle of an OWF. This has led to the development of various types of multi-objective decision-making models. In this section, the major types of these mathematical and spatial problems are categorized, and subsequent conclusions related to the spatial modeling processes are highlighted.

1.1. Key Aspects in Offshore Wind Energy Simulation Modelling

Choosing an optimal site is a complicated process based on various technical, environmental, and socio-economic parameters. Technical and geophysical factors must be examined, such as wind speed characteristics, wind power potential, water depths, sea substrate composition, and wave and current characteristics.

When considering the physical characteristics, a crucial step is estimating the wind power output and quantifying the associated uncertainty. According to different authors [6,7,8,9,10], the greatest uncertainty source in the energy output stems from selecting the optimal wind speed distribution. While most studies have derived wind power density assessments primarily from the Weibull distribution, multi-parameter distributions are generally more effective. These distributions tend to reduce errors in estimated wind power production, largely due to their ability to provide a better fit for extreme wind speeds (left and right tails of the distribution). Considering the annual energy production (AEP) calculation, methods like Weibull result in less accurate results and often underestimate the AEP by around 9% for a single year and around 4% for the entire operational lifespan of an OWF [6].

Understanding the previously mentioned characteristics of wind and power potential is essential for producing more accurate wind speed and power output forecasts or generating synthetic time series when long-term offshore wind in situ data is unavailable. To generate proper wind speed synthetic time series, it is necessary to select and consider the optimal marginal dependence of wind speed as also the diurnal and the seasonal wind speed behavior expressed by the long-term persistence, expressed as the average duration of wind speed in a given time interval for a certain location. Many models have been proposed in the literature, and the most applied ones are ARMA (Autoregressive Moving Average), Markov [11,12,13,14,15,16], and Wavelet models combined with Artificial Intelligence (AI) algorithms [17,18]. Although there are several methods for simulating arbitrary stochastic processes, they all have limitations and advantages in terms of the parameters used, the second-order dependence structure, and the wind’s intermittency preservation. Alternatively, a rigorous and general method was proposed by several studies [19,20,21] for the stochastic simulation of rain and wind fields, focusing on two key features. First, the climacogram stochastic tool (i.e., the variance of the average process vs. scale) for analyzing the variability of a natural process across different temporal scales, which can reveal the discretization, bias effects, and the inherent uncertainty via the Hurst parameter [20,22]. The second feature is the Symmetric-Moving-Average (SMA) stochastic simulation algorithm, which simulates a wind process while preserving marginal moments and significant stochastic properties, such as intermittency and long-term persistence [16,22,23,24,25], with the most known model in the literature to be the Generalized Hurst-Kolmogorov process.

1.2. Geographic Information Systems (GIS) in Offshore Wind Farm Planning

Emphasizing the spatial attributes, multiple studies have been conducted over the years based on geosciences, spatial modeling, and multi-criteria decision analysis for both cost and energy assessments. Indicative examples include the Netherlands, where a spatially explicit framework is established to assess the cost of OWFs along with the effects on the marine environment, the UK, China, Japan, India, South Korea, and the US [26,27,28,29,30,31,32], where GIS-based multi-criteria approaches were applied to map cost and energy-related indicators (foundation, transmission, installation, operation, and maintenance cost among with the annual energy production’s spatial distribution).

In addition, some robust and integrated spatial decision support systems (SDSS) have been developed by the northern European countries through the WINDSPEED program and the US [33,34,35,36] through a Web-Based Participatory GIS (PGIS) system in order to identify potential marine areas for OWF siting. In contrast, for the Mediterranean basin, significant research activity has been devoted to marine energy resources for Portugal, Italy, Egypt, and Greece [37,38,39,40,41,42], which are the first efforts to allocate space for OWF development using GIS and multi-criteria analysis tools.

One of the key characteristics of these studies is that they are based on monitoring and mapping the current or future state of potential marine areas. Nevertheless, several challenges remain unaddressed in modeling the spatial and temporal nature of energy systems at larger spatial scales, where it is essential to consider not only energy-related parameters but also geographic ones. According to References [43,44], this field of research has recently started to emerge, and developments in key fields, such as high-performance computing applications, big data analysis, and spatiotemporal modeling, allowed not only to improve model validation but also the development of geospatial tools with minimized uncertainties and spatial simplifications.

1.3. Spatial Optimization Models for OWF Site-Prospecting

Highlighting the challenges in geosciences, GIS played a pivotal role in spatial planning and site assessment procedures, although spatial optimization models for OWF planning are still in their infancy [45]. More holistic optimization approaches for OWF allocation and WT placement mostly focus on wind turbine level models and their layout configuration by examining the wind deficits (wake effects) [46,47] and inter-array cable cost minimization [48,49,50]. Furthermore, some challenging holistic projects developed over the last decade are based on multi-objective optimization for aerodynamics, the structural design of a WT, and the overall OWF cost minimization. The first holistic attempts applied with the Opti—OWECS, OWECOP, and OWFLO programs [51,52,53], having produced significant results in the field of wind turbine micro-siting problems and OWF configuration.

On the contrary, for site prospecting and suitability assessment processes, a few efforts are reported. For example, [45] used goal programming and decision-making methods and assessed pre-defined optional areas (9 proposed areas) for OWF siting. Their work accounted for the Life-Cycle Cost (LCC) analysis embedded in a multi-criteria approach. In addition, another study [54] developed an optimization process that produced promising results for identifying the most optimal and cost-efficient OWF locations in the UK. A comparison has been presented accordingly among the three state-of-the-art algorithms (NSGAII, NSGAIII, and SPEA2), but, similarly to [45], the selected approach was tested among pre-defined areas for future OWE deployment. In cases where no predefined sites are considered, a Genetic Algorithm (GA) can be employed [55] to predict candidate OWF locations based on criteria related to wind, depth, and distance; however, no consideration was given to shape or size configurations. Additionally, a raster-based Multi-Objective Turbine Allocation (MOTA) method was developed in [56] for onshore wind turbine site selection. Nevertheless, the focus of this approach was primarily on optimizing individual wind farm layouts rather than identifying optimal turbine clusters. At last, another study [57] used a K-Means clustering optimization approach based on random initialization of potential areas (point-based) upon different criteria (i.e., wind speed, water depth, and distance from shore) to identify favorable areas rather than OWF clusters for floating wind energy structures in Greece.

1.4. Critical Aspects for Spatial Optimization Site-Prospecting Procedures

To adopt a “tabula rasa” approach for site prospecting, it is essential to examine various factors and spatial considerations, including compactness and contiguity. Compactness is an indicator that belongs to the spatial optimization analysis that, in turn, applies diverse analytic and computational techniques involving geographical objects (cells, lines, points, areas-cells) [58,59,60,61,62]. Compactness refers to the shape and level of fragmentation among spatial units, therefore, it ensures continuity among spatial units is essential for achieving a compact area, meaning that neighboring geographical objects must be connected [63].

To address these issues, multiple research efforts have been applied in land use management and nature conservation, using either exact or heuristic spatial optimization techniques [62,64,65,66,67,68]. Numerous algorithms have been implemented to define and control compactness, but the most common and effective methods are based on (i) Integer Non-Linear Programming [69,70], (ii) Linear Integer Programming (LIP) neighbor methods [63,67], (iii) LIP using buffer zones cells [71], (iv) LIP using aggregated blocks [61,72], (v) minimization of shape index [60], and finally (vi) spatial autocorrelation (Moran’s method) [64].

Similarly, contiguity can be explicitly structured in spatial optimization models or implicitly accounted for in a solution algorithm [73,74]. Most explicit approaches are based on graph theory imposing network connectivity [68,70]. Furthermore, [56] highlights that incorporating the connectedness between spatial units within an optimization model often necessitates a substantial number of constraints and variables to guarantee parcels’ contiguity. Thus, graph-based models are commonly employed for this purpose. As a result, these approaches are typically confined to applications with limited spatial extents (i.e., a few dozen or hundreds of pixels) or are tailored to address specific land-use management problems.

1.5. Conclusions and Key Findings

As a concluding remark, several works [43,44,45] argue that GIS capabilities are not enough to solve some problems where the analytical and dynamic modeling aspects are crucial. Multi-objective modeling techniques for OWF planning and site prospecting are quite underdeveloped, and more integrated approaches must be taken into consideration. As a result, the present research aims to demonstrate an integrated spatial optimization scheme for the pre-development and planning stages of OWFs based on GIS structures and discrete optimization algorithms. The novelty of the proposed framework lies in the following:

• A spatial optimization framework is proposed as the first attempt for marine energy planning, considering some key spatial challenges and objectives, which are presented in detail elsewhere [42];
• An enhanced raster-based multi-objective optimization model is proposed (first applied in [62]), taking into account compactness by guaranteeing continuity via a Breadth-First Search cut algorithm;
• The climacogram stochastic tool is applied for the analysis of the long-term dependence structure with a focus on the variability of the wind power across a wide range of temporal scales [16,23]; and
• The applied Symmetric-Moving-Average (SMA) stochastic simulation algorithm with reduced parameters for computational efficiency, which is used for the simulation of the wind process while preserving, explicitly with the marginal moments, wind’s intermittency and long-term persistence as described in other studies in the literature (e.g., [16]).

The paper is organized as follows. In Section 2, the study area, data acquisition, calculation tools, and optimization techniques are introduced. The implementation and the experimental results are, respectively, given in Section 3, and finally, a summary, the conclusions, and the contribution of this work are reviewed in Section 4, along with suggestions and ideas for future research.

2. Materials and Methods

The outline of the proposed methodology is demonstrated in Figure 1 in a two-stage level of analysis incorporating the (i) offshore wind resource assessment and wind power stochastic simulation and (ii) multi-objective site-prospecting optimization algorithm. Central to this framework is a methodological implementation in an object-oriented programming language (Python) with the Gurobi optimization package (Gurobi Manual Version 2017 (Gurobi Optimizer Reference Manual Version 7.5)) for mixed-integer programming. Moreover, GIS software (ArcGIS 10.8.1) is used to store all available datasets in a geodatabase, to carry out all pre-processing procedures effectively, and to visualize the results. All calculation rules and sub-processes were automated for a spatial grid of a 5.5 km cell size (the spatial resolution of the UERRA reanalysis dataset). Therefore, all raster datasets are resampled to this spatial resolution and masked to a specific spatial extent.

2.1. Study Area and Data Management

In Figure 2, the study area is presented, located in the South–East Mediterranean Sea located between the Greek peninsula on the west and Asia Minor on the east, with coordinates 39° N 25° E. The North and Central Aegean Sea covers an area of 155,367 km². Greece extends to a little more than 130,000 km² and is characterized by its extended coastline of approximately 16,000 km long and by the existence of more than 3000 islands and islets, most of which are located in the Aegean Sea. According to other research works [42,75,76,77], the study area is characterized by large wind velocities (up to 9.4 m/s, 100 m above the mean sea level, mainly in the North-East and Central Aegean Sea). Also, unrestricted technical potential for offshore wind energy in Greece in 2030 based on average wind speed data is a bit less than 1000 TWh, highlighting the opportunities that arise for OWFs development in the Aegean Sea.

With its complex and highly irregular bathymetry profile, the Aegean Sea is one of the most physically challenging open water sites. Depths shallower than 250 m, considered technologically feasible for offshore wind energy, make up less than 50% of the total marine areas (see Figure 2).

Fixed-bottom structures like gravity-based, monopile, and jacket foundations are restricted to depths less than 70–80 m. However, floating structures can be placed at much greater depths (approx. 80–500 m.), although they have design challenges and technical constraints related to floaters, mooring lines, and anchors, which limit their depth limit [30,31,32]. The Northern Aegean is dominated by shallow shelf areas grading into a deep trench with three sub-basins: Sporades islands, Athos peninsula, and Lemnos Island. Deep basins, with a maximum depth reaching 1.5 km, depending on proximity to the West Anatolian fault, are near shallow areas. By contrast, the central and southern Aegean Sea have smoother bathymetric profiles. For instance, the Cyclades plateau reaches a depth of 400 m. The Central Aegean also includes deep depressions like the Skyros (east) and Chios Basins (west), as well as shallower regions on the west coasts of Lesvos Island, where depths exceed 250 m. The broader Cyclades and the southeastern Aegean Sea generally have shallower waters, but significant depth variations are seen near Ikaria and Samos islands [42]. Bathymetric data for the study were sourced from the Hellenic Navy Hydrographic Service (HNHS, https://hnhs.gr/en/, accessed on 12 November 2024) at a spatial resolution of 500 m.

The reanalysis wind data are retrieved from the UERRA (MESCAN-SURFEX) dataset of the ECMWF database (https://cds.climate.copernicus.eu, accessed on 12 November 2024), spanning 38 years (1982–2020). UERRA aims to recover sub-daily surface meteorological observations, which will be used as input and to assess future regional reanalysis products. The MESCAN-SURFEX data are available on a 5.5 km grid and at sub-daily timescales (6-h). The Numerical Weather Prediction (NWP) model used in regional reanalysis differs from those in global reanalysis. As a result, COSMO (Consortium for Small-scale Modelling) is applied in regional reanalysis, while the Hirlam Aladin Regional Mesoscale model (HARMONIE-ALADIN) is used to enhance spatial resolution further [78].

Numerous efforts of MESCAN-SURFEX reanalysis products are presented in studies in the literature [16,78,79], where high correlation measures indicated the precise capturing of time variability of the wind speed for the daily, monthly, and inter-annual scales. Although the data source comes with uncertainties since they are highly sensitive to bias issues, high-order moments and long-term behavior are effectively preserved. The statistical characteristics of the UERRA (MESCAN-SURFEX) dataset are illustrated in Figure 3.

2.2. Stochastic Simulation Model

The wind simulation process and the stochastic modeling scheme, based on the UERRA reanalysis records, include three key steps [16]. First, the climacogram stochastic tool [22] (variance of the average process vs. temporal scale) is employed to analyze long-term dependence, addressing variability across multiple timescales and handling discretization and bias effects in wind speed modeling. This is followed by applying the Symmetric-Moving-Average (SMA) stochastic simulation algorithm [80], which preserves all key statistical properties such as intermittency and long-term persistence through marginal moments [16,23]. The third step accounts for wind speed’s double periodicity (diurnal and seasonal) using an implicit scheme. When the periodicity is weak, it is suggested to divide periods, for example, by splitting a 12-month year into four seasons (December–February, March–May, June–August, and September–November) and a 24-h day into four 6-hour intervals (i.e., 00:00, 06:00, 12:00 and 18:00). For these parameters to be reliably estimated, the time-series must be long enough and of high quality, having discarded first the erroneous and missing values. Consequently, the stochastic generation process involves:

• Determining configurations and estimating parameters for each periodic cycle (diurnal and seasonal) and assessing persistence using the climacogram function.
• Applying generalized long-range dependence models for the wind simulation process, requiring the first four central moments, the Hurst parameter (H), the scale parameter (q) of the GHK model, and the synthetic time series length (i.e., N = 8760).
• Using the SMA scheme to preserve explicitly the first four central moments and the second-order dependence structure.

The GHK long-term dependence model is expressed via the climacogram, γ (m²/s²) (in Equation (1)) [16,25,80]. As a result, the simulated process (SMA) is expressed as the sum of the products of coefficients, and the long-term persistent behavior of wind is expressed by the standardized GHK model as

$$\mathit{\gamma }=\mathit{\lambda }/{(\mathbf{1}+\mathit{k}/\mathit{q})}^{\mathbf{2}-\mathbf{2}\mathit{H}}$$

(1)

$${\mathit{x}}_{\mathit{j}}={\sum }_{\mathit{j}=-\mathit{n}}^{\mathit{n}}{\mathit{a}}_{\left|\mathit{j}\right|}{\mathit{v}}_{\mathit{i}+\mathit{j}}$$

(2)

where λ is the standardized variance of the discretized stationary process, q is the shape parameter distinguishing the short-term from the long-term behavior, H is the Hurst parameter, ${\mathit{v}}_{\mathit{i}}$ is the white noise term, and ${\mathit{a}}_{\mathit{j}}$ can be analytically or numerically calculated up to lag n (see Algorithm 1).

Algorithm 1 Stochastic simulation model pseudo-code

Initialize the number of simulations N = 1000 and the total time series length L = 8760. Define the seasonal and diurnal cycles scale 1: Fit the BurrXII distribution of the Pareto-Burr-Feller (PBF) family to the entire reanalysis sample 2: Estimate for each periodic cycle the first m moments (16 m estimations for four seasons and four sub-daily 6-h intervals) 3: Estimate the m parameters of the PBF distribution for each cycle. 4: Transform each cycle through the inverse expression of the PBF distribution so that it follows the PBF distribution of the entire sample (Step 1). 5: Estimate the first four central moments via the PBF distribution. 6: Fit the GHK model with Hurst parameter H and scale parameter q, using the climacogram function: $\gamma =\lambda /{(1+k/q)}^{2-2H}$ 7: Estimate the coefficients ($a_j$) and the first four central moments of the white noise ($v_i$) of the SMA scheme. 8: Run the SMA model to generate the synthetic wind time series: $x_j$ = ${\sum }_{j=-n}^n{a}_{\left|j\right|}{v}_{i+j}$ 9: Transform back the synthetic wind time series using the parameters of the specific cycle (Step 2). Extract the total number of samples N with the synthetic offshore wind speed time series 10: Estimate the wind speeds at the height of 140 m (hub height) 11: Fit the wind speed time series (Step 10) to the power curve to estimate the mean AEP 12: Estimate the CV as a function of average wind speed, as extracted from the hindcasted UERRA Reanalysis data 13: Fit the GHK model, using the climacogram function, to the transformed wind power time series (N = 1000, Step 11) and estimate the upper and lower statistical limits of the variability. 14: Estimate the WPV as a function of the CV (Step 12) and the simulated power variability (Step 13)

The average marginal probability density function (PDF) for the standardized process is fitted and calculated based on the Burr XII distribution (Equations (3) and (4)), which includes the Weibull, gamma, and lognormal distributions [25]. It is also characterized as a distribution with a much heavier tail than the aforementioned and, thus, is found to be more appropriate for the wind process. The PBF distribution has two different asymptotic properties, i.e., the Weibull distribution for low wind speeds and the Pareto distribution for large ones [78], and with a PDF and cumulative distribution function (CDF), expressed as

$$f(x;b,c,a)={\displaystyle \frac{\mathit{b}\mathit{c}}{\mathit{a}}}({\mathit{x}/\mathit{a})}^{\mathit{b}-\mathbf{1}}[\mathbf{1}+({\mathit{x}/\mathit{a})}^{\mathit{b}}{]}^{-\mathit{c}-\mathbf{1}}$$

(3)

$$F(x;b,c,a)={\mathbf{1}-[\mathbf{1}+({\mathit{x}/\mathit{a})}^{\mathit{b}}]}^{-\mathit{c}}$$

(4)

with x ≥ 0, c and a as shape parameters, b as scale parameter, and for c*a ≥ r and c, a and b > 0. When c = 1, the PBF distribution becomes the Pareto II distribution, where a = 1, becomes a special case of the Fisk distribution, and for a → 0 becomes the Weibull distribution.

2.3. Annual Energy Production Calculation and Power Variability Quantification

The power curve for the Vestas V164-8.0 of 8 MW wind turbine (cut-in wind speed 4 m/s, rated wind power 8 MW at 13 m/s, cut-out wind speed 25 m/s) is used to estimate the Annual Energy Production (AEP) of an OWF. The mathematical relationship developed for Vestas V164-8.0 power curve [16] is presented in Equation (5). For each unit area, 25 wind turbines can be placed; therefore, wind turbine spacing is 7.5D, where D is the rotor diameter. With this interval, the influence of the wake is limited [42].

$$\mathrm{P}=\left\{\begin{array}{l} \mathbf{0} \mathit{V}\mathit{i}<\mathbf{4}\\ \mathbf{41.66} {\mathit{V}\mathit{i}}^{\mathbf{3}}-\mathbf{650} {\mathit{V}\mathit{i}}^{\mathbf{2}}+\mathbf{3858.3} \mathit{V}\mathit{i}-\mathbf{7600} \mathbf{4} \le \mathit{V}\mathit{i}\le \mathbf{7}\\ -\mathbf{25} {\mathit{V}\mathit{i}}^{\mathbf{3}}+\mathbf{775} {\mathit{V}\mathit{i}}^2-\mathbf{6500} \mathit{V}\mathit{i}+\mathbf{18,100} \mathbf{7}<\mathit{V}\mathit{i}\le \mathbf{11}\\ \mathbf{0.42} {\mathit{V}\mathit{i}}^{\mathbf{3}}+\mathbf{7063.1} \mathbf{11}<\mathit{V}\mathit{i}\le \mathbf{13}\\ \mathbf{8000} \mathbf{13}<\mathit{V}\mathit{i}\le \mathbf{25}\end{array}\right.$$

(5)

Wind power estimation has been carried out at a hub height of 140 m. Hence, the wind speeds need to be extrapolated to the required height for estimating the AEP per area unit (cell). For this purpose, the log-law (Equation (6)) has been used, where V(z) is the wind speed at a required measurement height z, V(${\mathit{Z}}_{\mathit{r}\mathit{e}\mathit{f}}$) represents the wind speed at the measured height ${\mathit{Z}}_{\mathit{r}\mathit{e}\mathit{f}}$, and ${\mathit{Z}}_{\mathbf{0}}$ denotes the roughness length. For simplicity reasons, constant sea surface roughness ${\mathit{Z}}_{\mathbf{0}}$ is assumed to be 0.002, based on previous studies [2,29], which is recommended for ocean surface and neutrally stable atmosphere.

$${\displaystyle \frac{\mathit{V}\left(\mathit{z}\right)}{\mathit{V}\left({\mathit{z}}_{\mathit{r}\mathit{e}\mathit{f}}\right)}}={\displaystyle \frac{\mathbf{l}\mathbf{n}(\mathit{z}/{\mathit{z}}_{\mathbf{0}})}{\mathbf{l}\mathbf{n}(\mathit{z}/{\mathit{z}}_{\mathit{r}\mathit{e}\mathit{f}})}}$$

(6)

Consequently, the annual energy production per unit marine area (MWh/km²/y) is estimated based on (Equation (7)) [26,29] considering an additional availability factor A of 90% and wake loss factor W of 10%:

$${\mathrm{A}\mathrm{E}\mathrm{P}}_{i,j}=(\sum _{t=1}^T\left({Pn}^{i,j,t}\right) A (1-W))nwt$$

(7)

where T is the total number of hours per year (8760), ${\mathit{P}}_{\mathit{n}}^{\mathit{i},\mathit{j},\mathit{t}}$ is the constant power output based on time t, and ${\mathit{n}}_{\mathit{w}\mathit{t}}$ is the total number of wind generators for each cell i, j.

Once the AEP is extracted for each simulation and cell, the climacogram function is applied to estimate the resulting wind power variability (WPV) per time scale in hours (i.e., 6, 12, 24, 48 … 8760). The resulting upper and lower bounds of the wind power variability delineate the resulting uncertainty of AEP (Equation (6)). To ensure that the resulting variability for each pixel includes the wind speed class characteristics (mean and variance), it is multiplied by the coefficient of variation (CV) [81,82] as expressed in Equation (8). Thus, to approximate the variability associated with the wind process, the trapezoid rule is applied as a method to calculate the value of the definite integral. The rule is based on approximating the value of the integral of a function g(x) by that of the linear function that passes through the points (a, g(a)) and (b, g(b)) expressed by

$$\mathit{WPV}=\left(\left(b-a\right){\displaystyle \frac{g\left(a\right)+g\left(b\right)}{2}}\right)CV$$

(8)

Notable is that these bounds represent the 1% (lower) and 99% (upper) confidence intervals of the estimated wind power output; however, alternate variability measures may be considered to quantify the estimated uncertainty.

2.4. Spatial Optimization Model

Formulating an optimization problem typically involves three fundamental steps: identifying decision variables, defining the objective functions, and establishing the problem’s constraints. To tackle multi-objective optimization problems, two primary approaches can be utilized: (i) the weighting method and (ii) the constrained method. The proposed algorithm employs a multiple parameters model using weighted summation, primarily for controlling spatial compactness, and is classified as an Integer Non-Linear Problem (INLP), solved using the Branch and Bound algorithm for superposition calculations.

For the integer programming formulation of the site-prospecting problem, a Boolean model is adopted where decision variables are linked to each potential OWF cell in binary form, indicating whether a cell is included in the solution (1) or not (0). This can be mathematically expressed as a combination of the maximization of the Annual Energy Production (AEP) and compactness index (CI) along with the minimization of the Wind Power Variability (WPV) and Depth Profile (DP) for each selected cluster of cells. To facilitate solving the problem, it was transformed into a multi-parametric optimization problem, aiming to minimize the final value of the objective function. This was achieved by inverting the values of the AEP index, making lower values desirable.

Algorithm 2 presents a pseudocode of the optimization process for identifying optimal areas of multi-megawatt projects (i.e., >500 MW) consisting of the predefined number of cells (N_CELLS) and clusters (N_CLUSTERS) distributed across the study area.

Algorithm 2 Spatial Optimization algorithm pseudo-code

Define all model parameters 1: Number of total clusters (N_CLUSTERS) 2: Number of total cells (N_CELLS) 3: Criteria weights $(W_k$) 4: Define the type of the optimization problem (i.e., minimization or maximization) 5: Normalization scheme applied for all criteria (i.e., 0–4 scale) Create model object 6: Add variables $X_{i,j}$ for any i in (0, rows + 1) and j in (0, columns + 1) 7: Add k “raster” surface values $C_{k,i,j}$ for any i in (1, rows) and j in (1, columns) 8: Add Contiguity Lazy constraint to check all optimal solutions 9: Contiguity check using BFS algorithm for every (i,j) in the current solution         If non-contiguous cells exist then: $\mathrm{ } \mathrm{sum}(X_{i,j}$) <= N_CELLS-1, for every (i,j) in the current solution 10: Add the Optimization Constraints $\mathrm{ } \mathrm{sum}(X_{i,j}$) = N_CELLS for any i in (0, rows + 1) and j in (0, columns + 1) $\mathrm{ } \mathrm{sum}(X_{0,j}$) = 0, for any j in (0, columns + 1) $\mathrm{ } \mathrm{sum}(X_{Rows+1,j}$) = 0, for j in (0, columns + 1) $\mathrm{ } \mathrm{sum}(X_{i,0}$) = 0, for any i in (0, rows+1) $\mathrm{ } \mathrm{sum}(X_{i,Columns+1}$) = 0, for any i in (0, rows + 1) 11: Set the Objective         If Minimization then: $\mathrm{ } X_{i,j}(W_1$ $C_{0,i,j}$ $+W_2$ $C_{1,i,j}$ $+W_3$ $C_{2,i,j}$ $+(W_4$ (4-sum(x[neighbors])))         Elif Maximization then: $\mathrm{ } {X_{i,j}(W}_1$ $C_{0,i,j}$ $+W_2$ $C_{1,i,j}$ $+W_3$ $C_{2,i,j}$ $+(W_4$ (sum(x[neighbors]))) 12: Create the solution array (in binary form, 0-No Solution pixels, 1-Solution pixels) 13: Save the solution array to raster(s)

2.4.1. Size, Data, and Complexity of the Study Area

The methodology uses raster files converted into arrays as inputs for the optimization algorithm, consisting of 50,588 discrete variables (binary), six linear constraints (two constraints linked with the total selected cells and four constraints for NoData cells handling), and one constraint (referred as lazy constraint) for the contiguity verification of the optimal solutions [42,62]. Each raster cell independently represents the (i) AEP, (ii) WPV, and (iii) DP values in array format (as extracted from the raster surfaces).

The compactness of a cluster is influenced by the surrounding neighborhood of each cell, with the objective of minimizing the total perimeter area (i.e., free edges) of the cluster. Specifically, the neighborhood of a cell is defined as a two-dimensional square lattice consisting of the central cell and its eight adjacent cells (Moore neighborhood) [62], as demonstrated in Figure 4. Among the total number of cells, NoData values within the raster surfaces are included but are managed by the model, as they do not contribute to the final optimal solutions. NoData cells can be borders, backgrounds, the sea, or other data considered to have invalid values. Consequently, the user can select or define constrained areas, mostly for planning or computational purposes. Figure 4 demonstrates how the algorithm handles the entire raster space, including the boundary cells. Considering an array Ci,j, an expansion of the total number of rows and columns is applied. Consequently, to be able to account for Moore’s neighborhood for the cells (e.g., above the coastline or an exclusion zone), two additional rows and columns are added, forming a final raster dataset Xi,j for i = i + 2 and j = j + 2 (blue array).

2.4.2. Decision Variables, Objectives and Constraints

Developing an optimization problem involves defining the decision variable(s), formulating objective function(s), and setting problem constraints (see the pseudo-code of Algorithm 2). These are selected based on the decision variables, expressed in raster-based formatted files (e.g., TIFF images) and transformed into arrays for input in the optimization algorithm. In this study, the following decision variables are selected:

• Maximum Annual Energy Production Index (AEP). The first objective is the maximization of the associated wind power output extracted from the stochastic wind field simulation in terms of the mean simulated values.
• Minimum Wind Power Variability Index (WPV). The second objective is based on minimizing the total variability index extracted from the wind power climacogram function and the coefficient of variation per cell.
• Minimum Depth profile index (DP). The third objective is based on minimizing the depths selected for potential OWF allocation. Therefore, different scenarios are considered by controlling the depth values.
• Minimum compactness index (CP). Compactness optimization is based on the normalized discrete compactness measure [83,84] consisting of a simplified approach applied to [64].

Considering the weighted criteria, the weighted sum method was employed as a fundamental and user-friendly approach for multi-objective optimization. As emphasized in the literature [85], when assigning weights for the a priori specification of preferences, it is crucial that each weight’s value is significant relative to other weights and their associated objective function. Consequently, only the relative importance of the objectives should be prioritized rather than the relative magnitudes of the function values. In this work, compactness and bathymetry are the most critical criteria for shape control and determining the technology used (i.e., fixed-bottom or floating foundations) in the optimal clusters. Accordingly, several additional runs were conducted with various weighting factors to identify the appropriate ranges for compactness (CP) and bathymetry (BP) weights, as well as the weights for the remaining criteria, AEP and WPV, and their inherent impact on the model configurations (see

Table A1. Optimization results in terms of the solutions’ quality and the computational performance for each scenario, with and without the Lazy Constraint activation.

Scenario	Cells	AEP	WPV	BP	CI	Lazy Constraints	Nodes Explored	Gap	Run Time	Obj. Value	Iterations	Contiguity
	(N)	Weight (0–1)				(Count)	(Count)	(%)	(s)	(Norm. cost)	(Count)	(Boolean)
1 (Lazy)	8	0.3	0.3	0.3	0.1	149/0	1	0	0.83	6.07	700	TRUE
2 (No Constraint)	8	0.3	0.3	0.3	0.1	0	1	0	0.68	6.01	474	TRUE
3 (Lazy)	8	0.2	0.2	0.2	0.4	2889/0	995	0	15	8.05	106,102	TRUE
4 (No Constraint)	8	0.2	0.2	0.2	0.4	0	1895	0	21.2	8.05	130,195	TRUE
5 (Lazy)	8	0.1	0.1	0.1	0.7	41,960/0	90,237	0	5678	10.02	17,234,568	TRUE
6 (No Constraint)	8	0.1	0.1	0.1	0.7	0	60,088	0	1,218	10.02	8,292,664	TRUE
7 (Lazy)	8	0.5	0.3	0.15	0.05	92/0	1	0	0.52	6.34	203	TRUE
8 (No Constraint)	8	0.5	0.3	0.15	0.05	0	1	0	0.46	6.34	169	TRUE
9 (Lazy)	8	0.5	0.3	0.19	0.01	155/0	1	0	Time limit	5.61	86	TRUE
10 (No Constraint)	8	0.5	0.3	0.19	0.01	0	1	0	0.0007	5.61	38	FALSE
11 (Lazy)	32	0.3	0.3	0.3	0.1	226,936/1057	108,785	0	346	23.73	245,565	TRUE
12 (No Constraint)	32	0.3	0.3	0.3	0.1	0	1	0	0.45	23.49	257	FALSE
13 (Lazy)	32	0.2	0.2	0.2	0.4	1214/0	267	0	8.75	23.25	20,534	TRUE
14 (No Constraint)	32	0.2	0.2	0.2	0.4	0	664	0	7.73	23.85	41,189	TRUE
15 (Lazy)	32	0.1	0.1	0.1	0.7	31,837/0	11,971	0	356	23.92	1,738,924	TRUE
16 (No Constraint)	32	0.1	0.1	0.1	0.7	0	12,923	0	224	23.92	1,389,533	TRUE
17 (Lazy)	32	0.5	0.3	0.15	0.05	471,234/3871	476,507	0	4287	25.98	2,876,765	TRUE
18 (No Constraint)	32	0.5	0.3	0.15	0.05	0	1	0	0.43	25.94	272	FALSE

in the Appendix A).

Figure 5 shows both alternative and optimal solutions for a continuous cluster of 5 cells within a hypothetical 5 × 6 raster area (30 cells). The four optimal alternatives, with a minimized total sum of −211, must be further evaluated based on compactness. Hence, the most compact cluster, determined by the smallest perimeter (calculated by free cell edges), has a perimeter of 12 for the first three scenarios (blue areas) and 10 for the fourth scenario. Continuity is guaranteed without applying contiguity constraints since the model minimizes fragmentation and prioritizes small perimeters. However, when compactness has a weight below 0.1, contiguity may not be ensured, leading to fragmented clusters (see Figure 5E). Additionally, increasing the compactness weight raises computational effort due to the non-linearity of the spatial optimization objective.

Contiguity control is based on a methodology that utilizes the Breadth-First Search (BFS) algorithm for traversing graph-based data structures [86]. This approach addresses the contiguity check as a graph-based problem, where each potential optimal solution is represented as a graph comprising edges and nodes (see Figure 6). Rather than being mathematically integrated into the optimization problem, contiguity is treated as a requisite property for any valid solution. According to a method following this approach [87], first is computing an optimal solution and then verifying whether all nodes of the graph within a contiguous solution have a connected path. If the solution meets all the constraints of the initial model, the optimization process is terminated; otherwise, the problem is resolved iteratively until an optimal solution is achieved. The proposed multiple factors optimization process uses weights for each criterion in the objective function and can be expressed as

$$Minimize=\sum _{(i,j∊C)}{x}_{ij}(w_{aep}{A}_{ij}+w_{wpv}{E}_{ij}+w_{bp}{B}_{ij}+w_v V_{ij})]$$

(9)

$$V_{\mathit{ij}}=4-(x_{i,j-1}+x_{i,j+1}+x_{i-1,j}+x_{i+1,j})$$

(10)

where C is the set of all raster ${cells}_{i,j}$ with total I rows and J columns, $w_V$ is the compactness weight, $w_{aep}$ and $w_{wpv}$ are the weights for the AEP and WPV respectively and $w_{bp}$ is the weight for the DP raster. For each ${cell}_{i,j}$, $A_{i,j}$ represents the value of AEP, $E_{i,j}$ the value of WPV, $B_{i,j}$ the DP and $V_{i,j}$ represent the total cell edges that contribute to the cluster perimeter; $X_{i,j}$ is a binary variable (0 or 1) to represent whether the ${cell}_{i,j}$ belongs to the final solution or not.

The optimization constraints in Equations (11)–(16) ensure that (i) the total number of selected optimal cells is equal to the desired size N (Equation (11)), (ii) none of the optimal cells exists at a pseudo-row or a pseudo-column of the raster (Equations (11)–(14)) and finally, (iii) the continuity (referred as Lazy constraint callback via the Gurobi optimizer) of the optimal solutions (Equation (16)). In particular,

$$\sum _{(i,j∊P)}{x}_{i,j} =N$$

(11)

where P is the set of all raster ${\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{s}}_{i,j}$ including two pseudo-rows and two pseudo-columns at the start and the end of the raster extent with a total I + 2 rows and J + 2 columns. For the neighborhood and the contiguity checks, the constraints are formulated as

$$\sum _{j=0}^{J+1}{x}_{0,j} =0$$

(12)

$$\sum _{j=0}^{J+1}{x}_{I+1,j} =0$$

(13)

$$\sum _{i=0}^{I+1}{x}_{i,0 } =0$$

(14)

$$\sum _{i=1}^{I+1}{x}_{i,J+1} =0$$

(15)

$$\sum _{(i, j ⫃\mathrm{S})}^{J+1}{x}_{i,j} =N-1$$

(16)

where S ⫃ C is an alternative solution with non-continuous raster cells (see the example in Figure 6).

3. Results

3.1. Stochastic Wind Process Simulation

The SMA-GHK model is known to successfully reproduce the marginal probability distribution, long-term persistence, and diurnal/seasonal periodicities (see [16] for detailed documentation).

Table 1. Simulated and estimated mean, variance, skewness and kurtosis between reanalysis and simulated data based on PBF distribution.

Marginal (Komotini)	Mean	Variance	Skewness	Kurtosis
Observed (UERRA)	3.679	5.942	1.281	5.751
Simulated	3.665	6.362	1.839	9.152
Marginal (Larissa)	Mean	Variance	Skewness	Kurtosis
Observed (UERRA)	3.713	7.971	1.714	8.511
Simulated	3.704	7.848	1.819	8.972
Marginal (Lemnos)	Mean	Variance	Skewness	Kurtosis
Observed (UERRA)	6.539	14.101	0.885	4.028
Simulated	6.506	13.916	0.961	4.322
Marginal (Cyclades)	Mean	Variance	Skewness	Kurtosis
Observed (UERRA)	7.271	12.211	0.606	3.361
Simulated	7.304	12.039	0.598	3.365

and Figure 7 and Figure 8 present the results of the synthetic wind speed and wind power time-series generation scheme regarding the marginal statistical properties of the UERRA and the simulated datasets (Table 1) and the second-order dependence structure (Figure 8), as a result of the Hurst parameter (H) and the slope (q). The simulated wind speed and wind power data consist of 1000 time series, each covering 1 year at a 1-h interval timestep.

The parameters related to the dependence structure via the climacogram (Equation (1)) are estimated for the entire study as λ = 0.18–78 (Figure 8a), q = 1 h, and the Hurst parameter (H) range in the interval 0.55–0.78 (Figure 8b), whereas for the marginal distribution (PBF) the shape (c,a) and scale (b) parameters are estimated in the range 1.2–2.5 and 1.6–50 for the shape and 2–20 for the scale parameter respectively. For example, when the shape parameter c increases, low wind speeds become more likely (left tail of the distribution) (Figure 7a). High shape parameter a values (Figure 7b) correspond to a heavy right tail, indicating extreme wind speeds, and higher scale b values (Figure 7c) represent greater wind speed values. Concerning the marginal characteristics, there is a strong agreement between the observed and simulated raw moments, including mean, variance, skewness, and kurtosis.

Table 2. Algorithm’s results for the running tests of 8 and 32 total optimal cells.

Scenario	No. of Cells	Compactness (CI) Weight	Nodes Explored	Simplex Iterations	Obj. Value	Run Time (s)
a	8	0.03	9286	26	2.94	62.5
b	8	0.1	9294	81	3.77	62.42
c	8	0.3	110	1978	5.71	4.04
d	32	0.03	5736	4509	16.42	63.96
e	32	0.1	5	400	17.69	2.46
f	32	0.3	158	982	19.51	3.82

provides a representative example for the regions of Komotini, Larissa, Lemnos, and the Cyclades (see Figure 8). Overall, the results indicate that the SMA-GHK model tends to overestimate the mean and slightly underestimate the variance in coastal areas while exhibiting the opposite trend in areas far from the shore. Although skewness and kurtosis show relatively minor variations (except for the Komotini area), their periodic cycles are effectively captured.

For the second-order dependence structure, the stochastic generation scheme supremely reproduces the observed long-term behavior of the wind power output (Figure 8). A historical wind speed time series is available, as extracted from UERRA dataset and simulated with N = 1000 synthetic series through the SMA-GHK process described in Section 2.1 and the statistical results presented in Table 1. The observed and simulated wind time series are then transformed to wind power output via Equations (4) and (5), and in Figure 8, the simulated and the real wind power variability and long-term persistence are identified for time scales larger than one year (i.e., 8760 h), and equals the half slope of the climacogram (Figure 8b), as the scale tends to infinity, plus 1. Analyzing the Hurst parameter (Figure 8a), the long-term persistence behavior of the wind power process is quantified and examined for the study area. Through the analysis of the entire study area, the Hurst parameter is estimated to be larger than 0.5 (0.52–0.65 for the Central Aegean Sea and 0.65–0.74 for the North Aegean Sea), indicating a long-term persistence that cannot be considered as white noise (i.e., H = 0.5). It is stressed that larger values of the Hurst parameter are associated with higher variability and uncertainty and, thus, less credibility and higher risk. Therefore, it is important to select the optimum areas for the wind installation turbines based not only on the average energy production but also on the lower associated variability.

Simulated wind power uncertainty from 1000 simulations is shown in Figure 8c. Decreased wind power variability at large timescales is noticeable in the wider offshore area of the Cyclades and some coastal parts (e.g., Larissa), where weak wind power potential also occurred. The upper and lower bounds of wind power variability (red lines), as shown in Figure 8c, outline the AEP uncertainty, and through Equation (8), the total WPV per pixel is extracted and incorporated into the optimization model.

3.2. Spatial Optimization Site-Prospecting

Initial tests focused on 8 and 32 cells with one cluster, using varying compactness weight values to evaluate the resulting areas’ shape, optimality scores, and computation time. All processing occurred on a desktop computer with an Intel Core i5 (2.3 GHz, 8 GB RAM) using Gurobi 7.5.2 to assess the proposed algorithm’s efficiency. The results and the computational effort are shown in Table 2 and Figure 9. For compactness control, the findings indicate that weight values below 0.1 lead to irregularly shaped optimal areas. Conversely, for weight values greater than 0.3, the shapes become more rectangular. Within the 0.1 to 0.3 range, irregularity decreases. Additionally, as the compactness weight drops below 0.1, the running time increases exponentially.

Figure 10 illustrates the results for the simulated wind speed, transformed to annual energy production (AEP) using Equations (4)–(6) and the wind power variability (WPV), as extracted from the coefficient of variation and the climacogram metric. All optimal areas were extracted based on varying numbers of cells, clusters, and bathymetry values (e.g., depths below or more than 60 m depth for bottom-fixed and floating foundations, respectively). The central and eastern regions of the study area, particularly the islands of the North Aegean Sea and the Cyclades, emerge as the most promising for offshore wind energy exploitation. In these areas, the estimated AEP ranges from approximately 650 to 900 GW annually, with the WPV values ranging between 0.5 and 2 (see Figure 10). In contrast, less favorable areas for wind energy production are located in the North-West Aegean Sea, including the Thermaikos Gulf, the Sporades Islands, and all coastal areas of the Attica region. Similarly, the eastern part of the Central Aegean Sea, particularly along the eastern coasts of Lesvos, Chios, and Ikaria Islands, also exhibits lower AEP values. These regions exhibit weaker offshore wind conditions, characterized as less suitable for large-scale offshore wind energy projects.

The optimization process, based on AEP and the WPV and the BP indices, was executed for three distinct scenarios to identify optimal areas for fixed-bottom and floating structures. The final results are depicted in Figure 10, illustrating three scenarios based on the number and size of clusters for OWF siting. The first scenario (Figure 10c) includes 3 clusters of 8 pixels, each with bottom-fixed foundations. The second scenario (Figure 10d) features 5 clusters of 8 pixels, each with floating structures. The third scenario (Figure 10e) comprises 2 clusters of 32 pixels each without any restrictions. To validate the model’s effectiveness in allocating continuous and compact patterns based on given weighting factors, different weights were tested for each objective. Thus, a weight of 0.4 was assigned to AEP and 0.2 to WPV, bathymetry, and compactness. The first two scenarios focused on areas with bathymetry below 70 m for fixed-bottom structures and between 70–200 m for floating structures, while the third scenario imposed no technical restrictions.

The most promising areas for bottom-fixed foundations are located in the North-East Aegean Sea, particularly around Lemnos Island and Alexandroupoli. Additionally, for floating structures and both floating and fixed-bottom structures (Scenarios 2 and 3), the broader central Aegean Sea shows optimal results, with annual energy production and wind power variability exceeding 900 GW annually. The spatial factors examined, such as site obstructions (e.g., depth restrictions) that limit the available area for OWE deployment, also indicate that water depth increases sharply, reducing the regions with shallow waters and significant wind power potential. However, suitable areas with water depths below 70 m were identified around the east coasts of Lemnos Island, in Alexandroupoli, the west coasts of Ikaria Island, and the east offshore areas of Central Evoia.

Although the average water depth for most European offshore wind farms is between 15–30 m and 10–20 km from the shore, this scenario is not fully applicable to the Aegean Sea. In many cases, such depths are within 10 km of the shore, creating proximity issues with potential conflicting onshore and marine activities, for example, near Alexandroupoli in the north and Lemnos and Ikaria islands in the northeast and Central Aegean Sea, respectively. This underscores the importance of careful spatial planning and stakeholder engagement to address competing uses and the inherent social impact.

4. Discussion

4.1. Stochastic Simulation and Uncertainty Quantification

The robustness of the SMA-GHK model is particularly evident in handling non-Markovian processes, such as wind speed, as also discussed in the literature [16,25]. By employing the climacogram and a multi-parameter distribution function (PDF), the model effectively captures the long-term persistence, the marginal moments, and the wind tails (i.e., the left tail of the distribution, which is associated with intermittency, and the right tail of the distribution that is associated with the extremes). These components are essential and must be accurately preserved across a broad range of temporal scales for OWE modeling.

The results for the simulated AEP are similar to those of other studies [16,75,76,77], indicating increased energy and wind power potential, particularly in the offshore areas extending from the North-East Aegean (around Lemnos Island) to the Central Aegean Sea (Cyclades Islands). This can be explained by the dry and strong Etesian winds (meltemia) in the summer period, specifically during August [16]. In addition to these findings, it is also reported that the observed annual and seasonal increases in wind power under various climate change scenarios (RCPs 4.5 and 8.5) for the Aegean Sea [88,89], all offshore areas along the western coasts of Lemnos, Lesvos, Chios, and Ikaria Islands show significant potential for future OWE exploitation. However, it should also be emphasized that the uncertainty associated with the estimated AEP presented in Figure 10a is not negligible due to the cumulative effects of various sources that cause underestimation. These include the UERRA wind speed inherent underestimation issues (see Table 6 in [16]), and the uncertainty introduced by the log-law extrapolation (approx. 1% per 10 m of height extrapolation [90]).

In addition to those findings, the climacogram was also chosen to explore the overall degree of unpredictability of the wind power potential by accounting for the Hurst parameter (H) and the coefficient of variation (CV). The use of the climacogram via the WPV measure is justified by the fact that CV has been characterized as a sensitive metric due to the presence of outliers or for distributions with long tails [82]. All examined processes display a Hurst parameter larger, on average, than 0.5, indicating a mild increased long-term uncertainty for specific areas in the North Aegean Sea (north, north-east offshore areas of Lemnos Island and all south coasts on northern Greece). In addition, the coefficient of variation (see Figure A1 in Appendix A) exhibited increased values in the same offshore areas as the Hurst parameter (H), similar to the results presented in previous studies [82,89] that used wind speed data from ERA-5 reanalysis and a multi-model ensemble of different Global Circulation Models (GCMs), respectively. These findings imply that only through detailed stochastic analysis may OWE deployment be more reliable and cost-efficient in all optimal areas.

4.2. Spatial Optimization Modelling Scheme

Analyzing the spatial distribution of the optimal clusters, a key advantage of the optimization scheme proposed in this study is its ability to generate solutions that closely align with users’ preferences, for example, the number of the optimal cells required, the compactness level (e.g., how rectangular is an optimal cluster) and the depth range (close or far from the shore) to be selected., The optimization model also employs a rigorous and systematic mathematical framework that avoids sub-optimal solutions (see the Gap Value of Table A1 in Appendix A), leading to optimal OWF sites in a reasonable amount of time. Moreover, users can modify the energy and depth raster surfaces by incorporating cost-related spatial information (e.g., foundation, transmission, and operation and maintenance costs), along with environmental, into the optimization process, similar to the criteria outlined in References [28,29,30,31,32,33,34,35,36,37,38,39,40,41]. This enables stakeholders and decision-makers to utilize the proposed spatial optimization methodology as a basis for investment decisions, integrating several key variables into the cost analysis of future OWE projects. Since no operational OWFs currently exist in Greece, model predictions cannot be validated against actual installations. Therefore, the accuracy of the model can only be assessed through comparison with suggested locations provided by the national authorities and previous research studies. Accordingly, the results presented in Figure 9 align with the potential organized development areas outlined in the draft National Programme for Offshore Wind Energy deployment in Greece (National Offshore Wind Farms Development Program) (see the map on pg. 120). These areas are located in the North-East Aegean Sea (Alexandroupoli), North Aegean Sea (western offshore areas of Chios Island), and Central Aegean Sea (western offshore areas of Ikaria Island). Similarly, the final optimal areas align with the regions where persistent patterns for optimal OWF placement are observed [40,41]. Specifically, the west coast of Ikaria Island is identified as an optimal location, consistent with Reference [40] (see Figure 6 in [40]). Also, the areas between Lesvos and Skyros Islands (Figure 10c,d—Clusters 1 and 2), as well as in Ikaria (Figure 10d—Cluster 5), also show alignment (see Figure 8 in [41]).

However, some limitations exist primarily in the model’s constraint setup. The spatial optimization results indicate that incorporating spatial contiguity as a constraint factor makes the model more challenging to solve using exact methods (in terms of the nodes explored and the number of iterations, Table A1 in Appendix A). Nonetheless, increased computational resources, broader availability of parallel computing, and advancements in solvers and methods to reduce complexity now enable the optimal solving of increasingly larger and more complex spatial problems. The results suggest that when a time limit is imposed, managing compactness through the objective function is the most efficient option to achieve optimal solutions with minimal computational effort. Conversely, if the shape of the optimal clusters is not considered and the optimality is based solely on cost objectives, contiguity control via the Lazy constraint becomes necessary, and this significantly increases the computational effort. While some problems were stopped after 15 min during the computational experiments, most solutions are designed to operate for extended durations, particularly when dealing with large-scale simulations that have low spatial resolution.

5. Conclusions

In this work, an integrated spatial optimization framework was described that represents the complex spatial relationships and interactions between multiple objectives for future OWF site-prospecting and marine energy planning. Based on this, some key tradeoffs inherent in OWF planning during the pre-feasibility stages were analyzed and discussed. The methodology can be easily utilized in other regions by applying the stochastic modeling scheme to extract wind power production and variability scenarios and the multi-objective constrained optimization algorithm for identifying optimal OWF clusters under different criteria.

In the first step, the wind time series generated with the SMA-GHK scheme were compared against UERRA reanalysis data based on (i) the marginal and periodic characteristics and (ii) the second-order dependence structure. The analysis of the stochastic wind simulation outputs suggests that, with the availability of extended periods of hindcasted reanalysis data, the proposed simulation scheme is capable of effectively preserving the marginal characteristics by applying the three-parameter PBF distribution and the GHK dependence structure fitted to historical wind time series via the climacogram estimator. Also, the utilization of the integrated WPV indicator, by employing the climacogram statistical tool, is emphasized for quantifying wind power uncertainty in comparison to the CV indicator. Thus, the robustness of the stochastic simulation scheme in Monte-Carlo experiments for sensitivity analysis is essential in renewable energy management to analyze predictability and optimal scenario planning under increased uncertainty and intermittency issues.

The multi-parameter constrained optimization algorithm demonstrated the proposed model’s effectiveness, efficiency, and potential in supporting pre-planning and decision-making processes for offshore wind energy (OWE) deployment. Hence, this framework offers planners and stakeholders a valuable tool for considering the inherent trade-offs between wind power output characteristics and cost-based spatial information, using straightforward implementation tools and objectives without resorting to complex spatial algorithms. Moreover, it is built on a rigorous and systematic mathematical formulation, ensuring optimal solutions by avoiding sub-optimal results and ensuring the identification of contiguous and compact offshore wind energy sites. Although initially developed for marine space management, this framework can be broadly applied as a powerful tool for location-allocation problems across various research fields, including the identification of future urban sprawl areas, extreme heat hot spots, and cumulative human impacts for Maritime Spatial Planning (MSP) processes.

Overall, this research presents valuable insights into spatial optimization for marine energy planning, yet several key findings pave the way for further exploration. These suggestions and challenges include (i) compactness control based on different shapes (i.e., lines, circles) for the resulted optimal areas, (ii) incorporating additional environmental (i.e., Marine Protected Areas—MPAs and essential marine and fish habitats), social (visibility impact) and cost-related criteria (e.g., wind turbine, foundation, installation, operation and maintenance cost and decommissioning costs) and test the spatial optimization scheme upon different cost indicators such as Levelized Cost of Energy (LCoE) or Net Present Value (NPV), as described in Reference [42], (iii) testing the stochastic simulation scheme on the Copernicus European Regional ReAnalysis (CERRA) dataset which showed increased accuracy in approximating in-situ measurements [91], and finally (iv) expand the proposed framework in a two-dimensional stochastic wind simulation framework. These directions not only aim to advance the robustness of the spatial optimization scheme but also address current limitations, laying a foundation for more comprehensive, site-prospecting and stochastic simulation analyses in related fields.