Mixing Renewable Energy with Pumped Hydropower Storage: Design Optimization under Uncertainty and Other Challenges

: Hybrid renewable energy systems, complemented by pumped hydropower storage, have become increasingly popular amidst the increase in renewable energy penetration. Such conﬁgura-tions are even more prosperous in remote regions that are typically not connected to the mainland power grid, where the energy independence challenge intensiﬁes. This research focuses on the design of such systems from the perspective of establishing an optimal mix of renewable sources that takes advantage of their complementarities and synergies, combined with the versatility of pumped hydropower storage. However, this design is subject to substantial complexities, due to the multiple objectives and constraints to fulﬁll, on the one hand, and the inherent uncertainties, on the other, which span over all the underlying processes, i


Introduction
In recent decades, the four pillars of the water-energy-food-land nexus have been significantly stressed.In an attempt to alleviate these stresses, systems that enable the development of different resources in the same area have been designed [1].Hybrid renewable energy systems (HRESs), which were first introduced in the 1970s, follow the same concept, since they combine multiple power sources over a limited geographical region, also sharing a common connection point.These have become an essential part of global energy production to address the limitations in terms of fuel flexibility, efficiency, reliability, emissions, and economics [2].HRESs generate electricity from two or more energy sources and are capable of switching between them when one is insufficient, thus reducing the inherited unpredictability of renewables.Moreover, they can also capitalize on the existing energy infrastructure and add components to help reduce costs, environmental impacts, and system disruptions [3].Overall, such systems can generate power based on the demand at any particular site, depending on the availability of resources, thus contributing to localized energy generation and supply management, and significantly reducing grid dependence [4].The optimal design of HRESs requires the most fitting, efficient, and reliable mix of energy technologies to meet users' needs.However, the asynchronous production of the intermittent renewable energy sources against the constant, yet highly fluctuating, energy demand is a major constraint against their large-scale integration in the energy grid [5].Interestingly, the absence of integrated large energy storage in the grid may result in a maximum of 80% wind penetration for a 100 kW grid, decreasing dramatically as the size of the electricity grid increases, diminishing to as low as 20% for a 10 MW grid [6].
For this reason, the rapid development of renewable energy systems and the associated technological advances (e.g., [7]) has necessitated their coupling with energy storage to balance the supply and demand.Notably, large energy storage systems have proven to mitigate wind power curtailment by 10% [8][9][10][11].Energy storage techniques are distinguished into four basic categories, according to their applications [12].Low-and medium-power applications (i.e., batteries, hydrogen fuel cells, and superconductors) concern isolated areas and are used in emergency terminals or individual electrical systems.On the other hand, network connections with peak leveling and power-quality-control applications (i.e., hydraulic systems, accumulators and flow batteries, and compressed air) are used in power systems that provide energy for larger regions.
In contrast to other renewables, hydroelectric energy, combined with water storage in reservoirs, enables the regulation of power production, energy recovery across water and oil transmission networks (e.g., [13,14]), and, even more importantly, energy storage [15].In fact, pumped hydro-storage (PHS) is the predominant and most reliable energy storage technology, accounting for more than 95% of the global cumulative energy storage capacity [16].Among its various benefits, Rehman et al. [17] highlighted several issues that PHS systems can address, which emerge from the high integration of wind power into the electricity network, such as (1) handling changes in network impedances due to wind farm connections to the grid and its effects on remote control signals, (2) the handling of harmonics created by the addition of wind to the grid, and (3) stability problems that may occur due to the dynamics behavior of wind farms connected to the grids.PHS systems are also not significantly influenced by the fluctuation in energy production occurring from the intermittent nature of renewables [18].In contrast, other well-established storage systems, namely, batteries, are prone to self-discharge effects when they remain in a loss power state for long periods of time, thus reducing their service life.Though PHS systems have been perceived as prohibitively costly, many studies have suggested their economic feasibility, especially in remote areas and islandic systems [19][20][21][22].
HRESs combined with a PHS system (which is in fact a water energy system) are mainly applicable in non-interconnected power systems, such as the ones found on islands.These are defined by geographical isolation, which constitutes their underwater interconnection to the main grid economically unviable, their low power demand, and high electricity production cost due to the importing of fossil fuels [23].In particular, Greece can greatly benefit from the aforementioned systems, since it consists of 29 non-interconnected islandic complexes [24], while its favorable weather conditions allow for high renewable energy potential.Moreover, the island's population increase during the summer leads to increased peak power demand, which further intensifies the need for energy independence.To date, only two HRESs have been installed on non-interconnected Greek islands (i.e., Ikaria and Tilos), as pilot applications, while there are also a few other systems under investigation, mainly in the research context (e.g., [25][26][27][28][29]).However, very few of these studies attempt to handle the major issue of uncertainty within the planning and design of such systems (e.g., [30]).
In general, PHS systems are configured as either open-or closed-loop systems.The first option involves the construction of an upper reservoir that is connected to a natural water source, which plays the role of a lower reservoir.On the other hand, closed-loop PHS systems refer to equally sized reservoirs without any external hydrological connections, recycling the same amount of water.When a PHS system is part of an HRES, the storage capacity of the reservoir system is dictated by the power tradeoff between the highly fluc-tuating and uncertain, as well, energy production by other renewables and the associated demand.Yet, conventional design practices follow much more simplified approaches in the context of PHS sizing.In particular, they aim at ensuring energy independence for a given time until the upper reservoir is fully discharged, during which weather conditions will not be favorable to sufficiently satisfy the power demand.However, this assumption is too conservative since it considers the worst-case scenario where energy is exclusively produced from the PHS system.Hence, it is expected to lead to reservoir oversizing, which in turn may threaten the technoeconomic feasibility of the overall project.
Recently, an HRES combining power production from renewables and energy storage through a PHS system utilizing seawater was proposed for the island of Sifnos [23].Its preliminary design followed the aforementioned conventional approach, resulting in a quite large reservoir capacity, if compared to the project scale.This case serves as a proof of concept to further expand a recently proposed stochastic simulation-optimization framework to the design of an HRES [31].
In this vein, this article aims at addressing a twofold objective.First, we investigate whether a more comprehensive optimization approach, driven by a detailed simulation model that accounts for the actual dynamics of water energy fluxes, can reduce the investment cost of the already proposed solution, by means of a smaller reservoir.Second, we explore the impacts of three key sources of uncertainty to the design procedure, involving two model inputs (wind velocity, as a physical process, and load demand, as an anthropogenic one), as well as the conversion of wind to electricity generation, through the power curves of the associated wind turbines.To our knowledge, there are very few attempts in the literature that aim to integrate and interpret all these different sources of uncertainty within the design of HRESs (cf.relevant review by [32]).A parallel objective, originating from the use of the sea as a lower reservoir, is to reveal the multiple technical challenges and resulting uncertainties, due to the use of seawater across multiple components of the PHS system (corrosion of the piping system, pumps and turbines, and risk of groundwater contamination in case of a leakage).

Study Area and Proposed Layout
Sifnos is a Greek island located in the Western Aegean Sea, Greece, in the Cycladic complex, with an area of 74 km 2 and a permanent population of 2755 residents, as per the 2021 census conducted by the Hellenic Statistical Authority.Sifnos attracts an average of 100,000 tourists during the summer months.Its energy needs are mainly covered by a 9.0 MW oil power plant, while renewables have a small share in the island's energy mix.Specifically, there is a 1.20 MW wind park and two photovoltaic parks of 0.20 MW cumulative installed power.According to an analysis of the island's energy profile for 2020, performed by the Hellenic Electricity Distribution Network Operator, the total energy demand was 17.3 GWh, while the hourly peak demand was 5.4 MW, occurring during the summer months.
A recent research study [23] investigated the development of an HRES in the northern part of the island (Figure 1) and its techno-economic perspectives.In particular, they represented the operation of an indicative layout, with different combinations of solar photovoltaic (PV) and wind turbine systems, combined with a PHS that utilized seawater.The upper reservoir was sited in a plateau, at an elevation of +320 m.This configuration was ideal, as there were no physical obstacles hindering the power production by the wind turbines and the PVs.Moreover, the favorable topography (i.e., steep, yet constant, terrain slope) facilitated the works concerning the conveyance system and the pumping station, thus requiring minimal technical interventions.Lastly, as social acceptance has a pivotal role in the implementation of renewable energy projects [33], we highlighted that the siting location was secluded, with the nearest settlement being 4 km away.A key outcome of the preliminary design analysis was the sizing of the PHS system, resulting in quite a large storage capacity of 1,100,000 m 3 .This capacity ensured energy autonomy for up to consecutive 16 days, starting from a fully charged state and without storing any excess energy during that time span.
Keeping the main layout of the aforementioned proposal, we attempted to revise some design quantities, and particularly the reservoir size, following a more integrated simulation-optimization procedure (Section 3).Furthermore, we embedded the key sources of uncertainty within the design procedure, which was formalized in stochastic terms (Section 4).

Configuration and Key Assumptions of the Simulation
This sub-section aims to specify the configuration of the HRES and outline the general assumptions of the underlying simulation model.The overall layout followed the one of [23], which was specified further after a preliminary analysis.Its components are illustrated in the schematic sketch of Figure 2, while its key characteristics are also summarized in Table 1.All essential design properties are specified, except for the reservoir's active depth and the number of PV modules, which are considered to be subject to the optimization.A key outcome of the preliminary design analysis was the sizing of the PHS system, resulting in quite a large storage capacity of 1,100,000 m 3 .This capacity ensured energy autonomy for up to consecutive 16 days, starting from a fully charged state and without storing any excess energy during that time span.
Keeping the main layout of the aforementioned proposal, we attempted to revise some design quantities, and particularly the reservoir size, following a more integrated simulation-optimization procedure (Section 3).Furthermore, we embedded the key sources of uncertainty within the design procedure, which was formalized in stochastic terms (Section 4).

Configuration and Key Assumptions of the Simulation
This sub-section aims to specify the configuration of the HRES and outline the general assumptions of the underlying simulation model.The overall layout followed the one of [23], which was specified further after a preliminary analysis.Its components are illustrated in the schematic sketch of Figure 2, while its key characteristics are also summarized in Table 1.All essential design properties are specified, except for the reservoir's active depth and the number of PV modules, which are considered to be subject to the optimization.A key outcome of the preliminary design analysis was the sizing of the PHS system, resulting in quite a large storage capacity of 1,100,000 m 3 .This capacity ensured energy autonomy for up to consecutive 16 days, starting from a fully charged state and without storing any excess energy during that time span.
Keeping the main layout of the aforementioned proposal, we attempted to revise some design quantities, and particularly the reservoir size, following a more integrated simulation-optimization procedure (Section 3).Furthermore, we embedded the key sources of uncertainty within the design procedure, which was formalized in stochastic terms (Section 4).

Configuration and Key Assumptions of the Simulation
This sub-section aims to specify the configuration of the HRES and outline the general assumptions of the underlying simulation model.The overall layout followed the one of [23], which was specified further after a preliminary analysis.Its components are illustrated in the schematic sketch of Figure 2, while its key characteristics are also summarized in Table 1.All essential design properties are specified, except for the reservoir's active depth and the number of PV modules, which are considered to be subject to the optimization.The wind turbines' tower height determined the wind speed values at the hub, which ere logarithmically distributed with respect to the elevation from the ground.Thus, two different wind turbines were chosen in order to minimize wind power curtailing due to cut-offs (either from very small or very large wind values).However, it was expected that the interaction between the large and small wind turbines (e.g., due to turbulence effects) would reduce the wind speed at the hub of the latter.The adjusted speed is calculated through the formula [34]: where V 0 is the freestream wind speed at the hub height level, a is an induction factor, k is a decay coefficient, L is the distance between the turbines, and D L is the large turbine's blade diameter.The a and k values are set equal to 0.10 and 0.038, respectively, as suggested in [34].In our case, we considered the use of four wind turbines, i.e., two large ones with 2.3 MW and two small ones with 0.9 MW of nominal power, respectively, distanced at 400 m.It is important to remark that the greater distance between the turbines might initially seem ideal, since the associated wind speed reduction is reduced; however, it does not necessarily result in optimal layouts [35].The wind power production for a given wind speed value is calculated by the analytical formula, introduced by Sakki et al. [36]: where P min is the minimum power produced during cut-in conditions, P max is the rated power, V wind is the actual wind speed, V min is the cut-in wind speed, V max is the cut-out wind speed, and a and b are the shape parameters.The two parameters were calibrated against the empirically derived power curve data provided by the manufacturer.For the wind turbine types used in our simulation, we obtained a = 2.25 and b = 20.The aforementioned formula can accurately describe the wind speed-to-power conversion process, as it considers turbine-specific characteristics, unlike the high polynomial order formulas that are commonly used to describe the wind power curve [37][38][39].
The following assumptions were also made for the PHS system: • The reservoir has a trapezoidal shape, and thus the storage and area curves are linear functions of elevation;

•
The intake is set to an elevation of 1.2 m from the upper reservoir's bottom to ensure sufficient capacity for deposit management; • The pump's power capacity is 6.0 MW and equal to the maximum potential surplus estimated as the difference between the total capacity of wind turbines (6.4 MW) and the minimum hourly demand (0.4 MW), occurring in winter during the night; • The turbine's power capacity is also 6.0 MW, which is slightly higher than the maximum hourly load (5.4 MW) in order to account for uncertainties, as discussed later;

•
The total efficiency values of the turbines and pumps are considered constant and equal to 0.85 and 0.80, respectively;

•
The penstock's length and diameter are 910 and 1.0 m, respectively, as specified in our preliminary design analysis.

Breakdown of the Simulation Model
The simulation model represents the system's operation in hourly time intervals, and is explained through a flowchart in Figure 3.At each time step, power accounting was performed by contrasting the power production by the two renewables (wind and solar) to the actual demand.If there were energy surpluses (P Net > 0), the PHS system was set to its charging phase, thus pumping water from the sea to the upper reservoir, provided that there was a sufficient storage capacity (S total < S max ).Similarly, if there were energy deficits (P Net < 0), the discharge phase was begun and the water was released downstream through the turbine, thus generating electrical energy, provided that there was available water stored in the upper reservoir (S total > V dead ).This process was repeated until the final step equaled the simulation length, which, in our case, was 20 years (which is the typical economic life of such projects).• The pump's power capacity is 6.0 MW and equal to the maximum potential surplus estimated as the difference between the total capacity of wind turbines (6.4 MW) and the minimum hourly demand (0.4 MW), occurring in winter during the night; • The turbine's power capacity is also 6.0 MW, which is slightly higher than the maximum hourly load (5.4 MW) in order to account for uncertainties, as discussed later;

•
The total efficiency values of the turbines and pumps are considered constant and equal to 0.85 and 0.80, respectively; • The penstock's length and diameter are 910 and 1.0 m, respectively, as specified in our preliminary design analysis.

Breakdown of the Simulation Model
The simulation model represents the system's operation in hourly time intervals, and is explained through a flowchart in Figure 3.At each time step, power accounting was performed by contrasting the power production by the two renewables (wind and solar) to the actual demand.If there were energy surpluses (PNet > 0), the PHS system was set to its charging phase, thus pumping water from the sea to the upper reservoir, provided that there was a sufficient storage capacity (Stotal < Smax).Similarly, if there were energy deficits (PNet < 0), the discharge phase was begun and the water was released downstream through the turbine, thus generating electrical energy, provided that there was available water stored in the upper reservoir (Stotal > Vdead).This process was repeated until the final step equaled the simulation length, which, in our case, was 20 years (which is the typical economic life of such projects).Notably, there was an ample amount of HRES analysis tools, whose functionalities were thoroughly explained in the comprehensive review of Sinha et al. [40].In our case, we formulated the entire computational procedure from scratch in the freeware RStudio programming environment (version 4.2.2).This allowed for embedding the specific peculiarities of the suggested layout and formalizing the simulation-optimization model in a stochastic (Monte Carlo) context.

Setup of the Optimization Problem
Techno-economic optimization is a crucial step within the development of complex systems, as it helps assess their feasibility and viability.This procedure involves systematically examining numerous technical, economic, and financial aspects, which can often be conflicting [41][42][43].
As already mentioned in Section 3.1, we considered as design variables of the system the active depth of the reservoir and the number of PV modules.The optimization of these variables was formulated as a multi-objective problem, by accounting for different performance criteria that were expressed in financial terms and weighted properly, in order to ensure the well-balanced operation of the hybrid renewable energy system.
Specifically, the overall performance measure sought to contrast the revenues from the electricity production with the costs of construction, purchase, installation, and maintenance.All the economic quantities were expressed in terms of the equivalent annual cost by considering a depreciation period equal to 20 years.The project implementation costs that were accounted for in the economic evaluation were associated with:

•
The purchase, installation, and maintenance of the electromechanical equipment (wind turbines, PVs, pumps, and turbines) and the conveyance system (GRP pipes); • Specific works associated with reservoir waterproofing.
A thorough cost breakdown for the HRES's construction and operation is presented in Table 2. Regarding the civil engineering works, their costs followed the typical pricing of respective projects in Greece.On the other hand, the cost of the pumps and turbines are calculated through the formula introduced by Aggidis et al. [44]: where C 0 = EUR 14,400, a = 0.56, β = −0.112,P is the nominal power and H is the gross head.By setting a power capacity of 6.0 MW (the same for the turbines and pumps) and a gross head of 320 m, we obtained a fixed total cost for the electromechanical components of the PHS system of about EUR 3,940,000.Finally, the maintenance costs of all individual works were set as 1% of the associated investment costs.A crucial requirement of electricity production systems is their reliability, which is a probabilistic concept, defined here as the frequency of load demand satisfaction.In our approach, we utilized the complementary definition of reliability, i.e., failure probability, which was quantified in economic terms and embedded in the optimization by means of an empirical penalty function.As illustrated in Figure 4, this function has a negative logarithmic shape and has been formulated after preliminary analyses, in order to ensure a desirable level of reliability.In this respect, the optimization aims at ensuring a financially attractive and technically sound investment.
Sustainability 2023, 15, x FOR PEER REVIEW 8 of 22 logarithmic shape and has been formulated after preliminary analyses, in order to ensure a desirable level of reliability.In this respect, the optimization aims at ensuring a financially attractive and technically sound investment.For the optimization, we utilized the evolutionary annealing-simplex algorithm (EAS) [45], which was a probabilistic heuristic global optimization technique that combined the robustness of simulated annealing in rough response surfaces with the efficiency of hill-climbing methods in convex areas.EAS is a well-established algorithm [46] that has been successfully utilized in a number of challenging optimization problems in the broader domain of water energy systems.

Results: Benchmark Scenario
The design optimization procedure results in an active depth of the reservoir equal to 2.88 m (total 4.08 m, by adding the freeboard of 1.20 m), and thus a storage capacity of 315,195 m 3 and a PV power capacity of 1.09 MW.Furthermore, the key metrics of the optimized benchmark scenario are presented in Table 3.We observed that the proposed solution ensured a quite satisfactory reliability level of approximately 95%; thus, the existing oil station will only have a complementary role in the island's energy mix by operating 5% of the time.We also underlined that the small capacity factor of the hydropower station (actual vs. potential energy production) did not indicate a reduced performance.In contrast, it revealed its pivotal supporting role in fulfilling the deficits by the other two renewables, especially during peak energy demand periods.As far as the other renewables' capacity factors were concerned, they were in line with the climatic regime of the study area.
Interestingly, the optimized reservoir size equaled less than a third of the one suggested in [23].However, the two solutions were not fully comparable since different assumptions were made in the overall modeling approach.This analysis will be used as the benchmark to formulate our stochastic optimization approach.For the optimization, we utilized the evolutionary annealing-simplex algorithm (EAS) [45], which was a probabilistic heuristic global optimization technique that combined the robustness of simulated annealing in rough response surfaces with the efficiency of hill-climbing methods in convex areas.EAS is a well-established algorithm [46] that has been successfully utilized in a number of challenging optimization problems in the broader domain of water energy systems.

Results: Benchmark Scenario
The design optimization procedure results in an active depth of the reservoir equal to 2.88 m (total 4.08 m, by adding the freeboard of 1.20 m), and thus a storage capacity of 315,195 m 3 and a PV power capacity of 1.09 MW.Furthermore, the key metrics of the optimized benchmark scenario are presented in Table 3.We observed that the proposed solution ensured a quite satisfactory reliability level of approximately 95%; thus, the existing oil station will only have a complementary role in the island's energy mix by operating 5% of the time.We also underlined that the small capacity factor of the hydropower station (actual vs. potential energy production) did not indicate a reduced performance.In contrast, it revealed its pivotal supporting role in fulfilling the deficits by the other two renewables, especially during peak energy demand periods.As far as the other renewables' capacity factors were concerned, they were in line with the climatic regime of the study area.Interestingly, the optimized reservoir size equaled less than a third of the one suggested in [23].However, the two solutions were not fully comparable since different assumptions were made in the overall modeling approach.This analysis will be used as the benchmark to formulate our stochastic optimization approach.

Issues of Uncertainty in Hybrid Renewable Energy Systems
Uncertainty has been a long-lasting issue in the planning, design, assessment, and real-time operation of HRESs, deriving from multiple drivers.Sakki et al. [31] divided such uncertainties into two main categories, i.e., exogenous (external) and endogenous (internal).The former category mainly refers to the inherent uncertainty of the system's drivers, whereas the latter refers to conversion processes and underlying modeling assumptions.In this vein, our study utilized the aforementioned stochastic optimization framework and applied it to the case of Sifnos by embedding three key sources of uncertainty (i.e., wind process, energy demand, and wind speed-to-power conversion).These uncertainties are briefly described in the following sub-sections, while the computational implementation is thoroughly explained in Section 5.

Wind Process Uncertainty
Hydrometeorological processes are considered one of the main exogenous uncertainties of HRESs due to the intermittent nature of renewables.Specifically, these concern wind processes (i.e., wind velocity and wind direction) that involve Aeolic systems, hydrological processes (i.e., precipitation, streamflow, etc.) that are associated with hydropower systems of all scales (from large reservoirs to small run-off-river plants), as well as solar-related processes (i.e., solar radiation, cloud cover, etc.) that are associated with photovoltaic energy.
Regarding wind power, which was one of the key components of our uncertainty analysis, the mainstream approach to represent the variability of wind velocity was to generate synthetic data at coarse time scales (e.g., monthly or daily) by applying theoretical distribution functions following the statistical regime of the historical data.However, it is well known that wind speed exhibits significant fluctuations, even across very small scales, thus requiring a much finer modeling resolution.Furthermore, simple statistical tools fail to capture a key aspect of all hydrometeorological processes, namely, the dependencies on time and space.In this respect, the most appropriate approach for representing the full regime of such processes (including wind) are stochastic models, applied at fine scales (typically hourly).
In this context, several researchers performed comprehensive analyses of hydrometeorological data worldwide and proposed suitable stochastic methods for simulating the input drivers of renewable energy systems.For instance, Palma et al. [47] presented a novel methodology to facilitate the selection of a proper time-series generation model for renewable energy sources, providing a set of indicators to verify the selected model's accuracy.Furthermore, other studies also emphasized another key property of such processes, which was referred to as long-term persistence and was associated with the changing climate [48-50].

Energy Demand Uncertainty
Energy demand is an even more complex source of exogenous uncertainty, as it does not only depend on physical climatic processes (mainly temperature), but also on highly unpredictable environmental and socioeconomic factors [51].Uncertainty in electricity demands, both for simulation and forecasting purposes, is a topic of high interest in the literature.For instance, a study conducted by Cabeza et al. [52] showed that past energy projections for the energy demand of OECD countries were either consistently overestimated or underestimated.This phenomenon intensified even more on the islands due to the intensive seasonal power demand variations.Thus, it is essential to accurately represent energy demand in probabilistic terms while designing HRESs.Warren et al. [53] also highlighted three main constituents of uncertainty in energy demands: (1) inherent randomness in the way electricity is consumed, (2) modeling and estimation errors, and (3) uncertainty in the model inputs.Finally, Islam et al. [54] presented various models for short-, medium-, and long-term energy demand modeling practices under uncertainties and metrics to effectively evaluate their accuracies.

Wind-to-Power Conversion Uncertainty
As the penetration of wind power in energy systems increases, several concerns about the uncertainty in wind power generation have been raised.Uncertainty in wind power system operations can be categorized between discrete and continuous disturbances [55].The discrete disturbances causing equipment failure, involving generators and transmission lines, were accurately demonstrated in an analysis performed on wind turbine generators by Rezamand et al. [56].The results indicate that the reliability of wind turbine generators (WTGs) can decline to as low as 67.9% after seven years of operation.The continuous disturbances, which include parameters of the unit commitment problem that vary smoothly (e.g., electricity demand and wind power production), were described in the previous sub-section.
Another type of uncertainty to be considered, which is also emphasized in our study, involves the wind turbine power curves (WTPCs), which have a pivotal role in the context of wind power simulation and forecasting, wind turbine condition monitoring, and in the estimation of wind energy potential [57].Given that wind turbines are commercial products, WTPCs are provided by the manufacturers as a graph or a set of points with a typical wind speed discretization of 0.5 m/s.In fact, these are theoretical relationships exhibiting a wind turbine's standard and experimental behavior [58].However, a wind turbine operates in complex and variable conditions, which deviate significantly from the stable experimental conditions under which manufacturers test them.Thus, the herein referred to as "theoretical" WTPCs cannot accurately represent the actual behavior of wind turbines that operate in the field [59].The deviation of wind turbine on-site behavior from the theoretical power curve was thoroughly analyzed in the work of Antoniou et al. [60], whilst several real-world examples were also presented in [61].For this reason, an ample number of deterministic and probabilistic models were developed to produce WTPCs that resembled actual operating conditions.Recently, the focus has shifted towards the latter since deterministic models provide fixed relationships between wind speed and power generation, failing to reveal the variating and dynamic power generation process [62].A novel probabilistic WTPC model worth mentioning is the one developed by Yan et al. [63], which considers various model inputs (pitch angle and wind direction) based on three non-parametric algorithms.
A major issue of uncertainty involves the operation of wind turbines in a high-windspeed region.According to the theoretical model, the turbines are forced to interrupt their operation in the so-called cut-off wind speed of 25 m/s to prevent damage due to extreme mechanical loads.The turbines start operating only after the average wind speed reaches a value lower than the cut-off speed, a process which is also referred to as hysteresis [64].However, due to the stochastic nature of wind and the uncertainty of the theoretical WTPCs, it is not realistic to consider the threshold of 25 m/s as a strict rule for estimating the downtime of turbines [65].Thus, to prevent frequent shutdowns and restarts, the soft cut-out strategy was implemented by extending the maximum admissible wind speed up to 30-32 m/s, without an abrupt shutdown, through controlling the pitch and generator in order to gradually decrease energy production [66].Multiple studies have addressed the optimal control of wind turbines during high wind speeds to minimize the uncertainties derived from the wind.Jelavic et al. [67] produced a soft cut-out strategy worst-case scenario algorithm that did not significantly increase the fatigue loads.Astolfi et al. [68] performed a SCADA data analysis, extending the power curve of a wind turbine farm in the high-speed region, and concluded that the simulated energy improvement was 0.62%, which equaled 1.80% of the wind farm's total production.Lastly, Castellani et al. [66] extracted operational data from a 2.3 MW wind turbine, which was then simulated to work with the soft cut-out strategy, producing 1.02 MW more energy than its initial operating state.
A final source of uncertainty with respect to actual wind power generation, in contrast to the theoretical behavior of the power curve, was associated with monitoring errors.The output power observations were acquired through Supervisory Control and Data Acquisition (SCADA) systems.However, such systems often contain abnormal data and outliers, occurring from wind curtailment, maintenance, or other uncontrollable factors [69].Consequently, these data need to be cleared by appropriate classification algorithms in order to eliminate their adverse effects and improve wind power predictability [70].

Incorporating Uncertainty in the Simulation
Following the uncertainty aware framework by Sakki et al. [31], we ran the design optimization problem in the Monte Carlo context.This involved the formulation of many scenarios (specifically, one hundred), which allowed for providing equally probable sets of optimized solutions instead of a unique one.
In Section 1, we outlined the three key sources of uncertainty that were embedded in the aforementioned Monte Carlo simulation-optimization approach, categorizing them into two external processes (wind velocity and energy demand) and an internal one (wind power conversion to electricity).The first two sources of uncertainty ere integrated in the simulation by means of synthetic time series that were generated through the anySim package [50], which provided a suite of stochastic models for the simulation of both stationary and cyclostationary processes (in univariate or multivariate modes) that may follow a wide spectrum of distributions.This allowed for generating synthetic time series with the desired marginal and stochastic properties, as reflected in the associated historical data [71][72][73].In our work, we utilized this package, which ran in the R programming environment, to generate 100 synthetic time series of hourly wind velocity and hourly energy demand (one for each scenario), for a 20-year horizon, which is the typical economic lifespan of an HRES.Their probabilistic and dependence regime reproduced the one of the corresponding historical data, the length of which (ten years) was rather short for a proper representation of their actual variability and, thus, uncertainty.
Regarding the uncertainty in the wind velocity to power conversion, as mentioned in Section 4, wind turbines operate in complex and variable conditions, thus the power curves provided by the manufacturers did not accurately represent their on-site behavior.In this vein, we followed a Monte Carlo simulation approach, considering that the shape parameters a and b of the empirical Equation (2) were random variables.Specifically, for both turbines, the output power for a given wind velocity value was sampled from a normal distribution with a~N(2.25,0.0225) and b~N(20, 0.016).This distribution and its associated parameters (i.e., mean and standard deviation values) were appropriately selected after some trials, so that the average deviation from the theoretical curve was not greater than 15% (Figure 5).This deviation was within reasonable margins, since Lira et al. [74] stated that the power curve uncertainty was approximately 10%.
Furthermore, we also implemented the soft cut-out strategy, allowing the wind turbine to operate in the high-wind-speed region (Figure 5).Therefore, we considered a linear reduction in the turbine's power for wind speed values between 25 and 30 m/s.The slope, z, of this linear equation was also represented as a normal process, with z~N(150, 225) for the large turbines and z~N (45,100) for the small ones.Furthermore, we also implemented the soft cut-out strategy, allowing the wind turbine to operate in the high-wind-speed region (Figure 5).Therefore, we considered a linear reduction in the turbine's power for wind speed values between 25 and 30 m/s.The slope, z, of this linear equation was also represented as a normal process, with z~N(150, 225) for the large turbines and z~N (45,100) for the small ones.

Results of Monte Carlo Scenarios
Table 4 summarizes the key statistical characteristics of the basic outputs of the 100 optimization scenarios, which are expressed in terms of mean, standard deviation, and three typical quantiles, i.e., 10, 50, and 90%.We remembered that the two design variables to be optimized were the reservoir active depth and total power capacity of the PVs.The conventional approach of Section 3 resulted in the optimized values of 2.88 m and 1.09 MW, respectively.However, under the Monte Carlo framework, the reservoir depth and, consequently, its storage capacity, exhibited significant variability, which was propagated in the key performance metrics of the system (net profit and reliability).In contrast, the variability of the PV capacity was negligible.In this context, the range of uncertainty with respect to renewable energy production mainly reflected the stochasticity of the wind process and its conversion through the probabilistic power curves of the two turbines.We observed that while the installed wind power capacity remained constant (i.e., 6.4 MW), the capacity factors of the two machines varied (although, not significantly) as a result of this stochasticity.On the other hand, since for all scenarios we applied the same set of solar radiation data, the capacity factor of the PVs remained constant and equaled 0.21.

Results of Monte Carlo Scenarios
Table 4 summarizes the key statistical characteristics of the basic outputs of the 100 optimization scenarios, which are expressed in terms of mean, standard deviation, and three typical quantiles, i.e., 10, 50, and 90%.We remembered that the two design variables to be optimized were the reservoir active depth and total power capacity of the PVs.The conventional approach of Section 3 resulted in the optimized values of 2.88 m and 1.09 MW, respectively.However, under the Monte Carlo framework, the reservoir depth and, consequently, its storage capacity, exhibited significant variability, which was propagated in the key performance metrics of the system (net profit and reliability).In contrast, the variability of the PV capacity was negligible.In this context, the range of uncertainty with respect to renewable energy production mainly reflected the stochasticity of the wind process and its conversion through the probabilistic power curves of the two turbines.We observed that while the installed wind power capacity remained constant (i.e., 6.4 MW), the capacity factors of the two machines varied (although, not significantly) as a result of this stochasticity.On the other hand, since for all scenarios we applied the same set of solar radiation data, the capacity factor of the PVs remained constant and equaled 0.21.

Insight into the Trade-Off between Reservoir Size and Overall System Profit
The analysis of the model outcomes under uncertainty revealed the existence of a trade-off between the key design quantity, i.e., the active depth of the reservoir, which dictated its size, and the overall profit of the water energy system.Since this trade-off embedded all kinds of uncertainties that were included in the simulation-optimization procedure, its quantification and interpretation presupposed the use of advanced stochastic and statistical tools.This objective was implemented via a two-step approach.
First, we investigated the marginal behavior of the two random variables by assigning suitable statistical distributions, specifically normal and log-normal ones, to the reservoir active depth and the mean annual profit, respectively (Figure 6).In this context, we utilizes the fitdistrplus R package [75], which provided functions for fitting a wide range of univariate distributions to different types of data (continuous censored or non-censored, and discrete as well), also supporting different parameter estimation procedures.In our case, for both variables, we applied the moment-matching estimation method, which involved finding the values of the model parameters that made the data's sample moments equal to the model's corresponding population moments.
utilizes the fitdistrplus R package [75], which provided functions for fitting a wide range of univariate distributions to different types of data (continuous censored or non-censored, and discrete as well), also supporting different parameter estimation procedures.In our case, for both variables, we applied the moment-matching estimation method, which involved finding the values of the model parameters that made the data's sample moments equal to the model's corresponding population moments.
After selecting the marginal distributions for the two variables, we investigated their dependencies by fitting a Gaussian copula, as shown in Figure 7. Notably, the two variables exhibited a quite strong negative correlation, as indicated by the associated Pearson's coefficient (R), which equaled −0.75.The copula approach enabled the formulation of multivariate non-normal distributions by combining given non-normal marginal models, only with dependence patterns.Gaussian copulas are very flexible and have been broadly used for the modeling of dependent variables of any type [76,77], as well as for predictive uncertainty assessment studies.
In our case, the copula-based approach provided greater insight into the trade-off between the optimized reservoir depth and the associated profit under uncertainty.In particular, for a given depth, and thus a known reservoir size, which is an engineering decision, we may specify a probabilistic range of the anticipated annual profit of the system.On the contrary, a desirable profit threshold corresponds to a range of potential design options.In this vein, a conservative design approach would favor the implementation of a large reservoir (upper quantile) and vice versa.Conclusively, the copula graph can be utilized as a decision support tool for both engineers and from the perspectives of investors and stakeholders.After selecting the marginal distributions for the two variables, we investigated their dependencies by fitting a Gaussian copula, as shown in Figure 7. Notably, the two variables exhibited a quite strong negative correlation, as indicated by the associated Pearson's coefficient (R), which equaled −0.75.The copula approach enabled the formulation of multivariate non-normal distributions by combining given non-normal marginal models, only with dependence patterns.Gaussian copulas are very flexible and have been broadly used for the modeling of dependent variables of any type [76,77], as well as for predictive uncertainty assessment studies.

The Challenge of Seawater
One of the major obstacles for the development of HRESs in small, non-interconnected islands, which take advantage of PHS systems, is the scarcity of water resources, which is a long-lasting issue, especially during summer.A potential solution can be the deployment of PHS systems that utilize the sea as the lower reservoir, as proposed for the case of Sifnos island.The Okinawa Yanbaru Seawater Pumped Storage Power Station [78] was the first experimental pumped storage facility in the world that used seawater for energy storage.Its construction lasted from 1987 to 1999; yet, it was dismantled in 2016.This was because Okinawa's projections overestimated the load demand growth, constituting the plant's operation as no longer profitable.
Over the years, there have been several studies suggesting the use of PHS systems that utilize the sea as a lower reservoir in islands, including Hawaii, Ireland, and the Azores [79][80][81].In a more recent study, Ali and Jang [82] examined the case of Deokjeokdo, South Korea, and found that an HRES coupled with a seawater PHS system, apart from providing energy independence, also led to a cheaper initial investment cost than using batteries as energy storage.However, the inclusion of a seawater pumped storage system in the HRES introduced various technical challenges, mainly related to corrosion.In our case, the copula-based approach provided greater insight into the trade-off between the optimized reservoir depth and the associated profit under uncertainty.In particular, for a given depth, and thus a known reservoir size, which is an engineering decision, we may specify a probabilistic range of the anticipated annual profit of the system.On the contrary, a desirable profit threshold corresponds to a range of potential design options.In this vein, a conservative design approach would favor the implementation of a large reservoir (upper quantile) and vice versa.Conclusively, the copula graph can be utilized as a decision support tool for both engineers and from the perspectives of investors and stakeholders.

The Challenge of Seawater
One of the major obstacles for the development of HRESs in small, non-interconnected islands, which take advantage of PHS systems, is the scarcity of water resources, which is a long-lasting issue, especially during summer.A potential solution can be the deployment of PHS systems that utilize the sea as the lower reservoir, as proposed for the case of Sifnos island.The Okinawa Yanbaru Seawater Pumped Storage Power Station [78] was the first experimental pumped storage facility in the world that used seawater for energy storage.Its construction lasted from 1987 to 1999; yet, it was dismantled in 2016.This was because Okinawa's projections overestimated the load demand growth, constituting the plant's operation as no longer profitable.
Over the years, there have been several studies suggesting the use of PHS systems that utilize the sea as a lower reservoir in islands, including Hawaii, Ireland, and the Azores [79][80][81].In a more recent study, Ali and Jang [82] examined the case of Deokjeok-do, South Korea, and found that an HRES coupled with a seawater PHS system, apart from providing energy independence, also led to a cheaper initial investment cost than using batteries as energy storage.However, the inclusion of a seawater pumped storage system in the HRES introduced various technical challenges, mainly related to corrosion.This section aims to address the most crucial ones, ensuring the proposed system's technical feasibility.

Conveyance System
Corrosion prevention can vary depending on whether the components are subject to high or low flows.This protection categorization was implemented in the Okinawa project [83].The conveyance system of the latter utilized pipes fabricated of fiber-reinforced plastic (FRP) with rubber joint seals, which are resistant to both seawater corrosion and high pressure.Generally, glass-reinforced polyester (GRP) was chosen for smaller pipe diameters [84].GRP is a composite material that consists of a polymer matrix and glass fibers.The former is usually an epoxy, vinyl ester, or polyester thermosetting resin, acting as a binder for the fibers, while the latter adds strength to the composite [85].Its chemical structure is inert to seawater and its smooth surface reduces hydraulic losses, exhibiting a quite satisfactory friction factor of 0.030.However, this material's internal water pressure decreases as the pipe's diameter increases.Nevertheless, GRP pipes can be combined with steel pipes, if higher pressures are expected [84].
Another important issue regarding the conveyance system is the unwanted accumulation of marine organisms immersed in the sea, which is also referred to as marine fouling [86].In the past, various methods were used to prevent that phenomenon, such as chlorine production by an electrolysis unit and its disposal in the penstock [87].However, it was found that such practices not only threatened the surrounding ecosystem, but also deposited byproducts (i.e., manganese) [88].While less harmful methods that do not produce biocidal intermediates have been investigated, i.e., electrochemical degradation [89], natural biocide-based, and non-stick coatings [90,91], further research is required in order to provide solutions that ensure both the surrounding ecosystems' sustainability and the project's techno-economic feasibility.

Electromechanical Equipment
Electromechanical equipment refers to the turbines and pumps used in the PHS system.In HRESs, the selection of pumps is a particularly important procedure, as there is a wide variety available on the market, serving various applications.For a project of this scale, centrifugal pumps are recommended, given that they have been previously used for marine applications [92].Their basic limitation is associated with the intermittent production of renewables, meaning that when weather conditions are not favorable, it may be impossible to develop sufficient heads to fill the reservoir.This issue was addressed by Manfrida and Secchi [92], who suggested a configuration, including the installation of pumps in parallel, with different flow rates.
Similar to the corrosion of pipes, the erosion of electromechanical equipment poses another technical challenge.Francis and Hebdon [93] distinguished several types of corrosion from seawater that could affect stainless steel (SS) electromechanical equipment into the following ways:

•
Crevice corrosion, which is the most ordinary form of corrosion, is initiated by changes in the local chemistry within a crevice.It is usually associated with a stagnant solution in the micro-environments that tends to occur in crevices.In seawater pumps, crevices can be found where seals and impellers are fastened to the shaft and flange faces are cast in for pipework connections; • Erosion corrosion can occur from the seawater's rapid flow rate; • Cavitation occurs when a fluid's operational pressure drops below its vapor pressure and causes gas pockets and bubbles to form and collapse.This common phenomenon occurs when a pump operates outside its normal design parameters.The formed bubbles erode the steel; • Corrosion fatigue derives from the combination of alternating or cycling stresses in a corrosive environment, mainly affecting seawater pump shafts.
In general, alloys of stainless steel (i.e., chromium, nickel, and molybdenum) have a high pitting-resistance equivalent number (PREN), which indicates its corrosion resistance.According to the Norwegian Standards, a PREN value of around 40 provides sufficient corrosion resistance to seawater [94].

Groundwater Degradation Due to Seawater Effects
The last challenge introduced by seawater concerns the possibility of ground contamination by salt water (either by leakage or wind transport).Specifically, groundwater salinization is a crucial environmental issue, especially in areas with karstified underlying geological formations (this was also one of the major shortcomings of the proposed design for Sifnos).Thus, the waterproofing of the reservoir is considered essential while dealing with seawater.In this context, high-density polyethylene geomembranes (HDPEs), which are resistant to chemicals and ultraviolet radiation, are a potential solution, combined with a drainage system that detects possible leakages [78].On the other hand, the construction of an embankment around the reservoir can prevent seawater transport via the wind [84].

Conclusions
In the present study, we developed a simulation-optimization procedure for the design of HRESs in a stochastic setting and assessed it within the preliminary analysis of a proposed scheme for a non-interconnected Aegean island (Sifnos).The HRES model was based on an hourly time-scale simulation approach, accounting for a detailed representation of all important meteorological, energy, and water fluxes.With regard to this, the optimization method opted for a robust and techno-economically sustainable outcome, thus maximizing the project's mean annual profit until its full depreciation and simultaneously ensuring high reliability in the load demand satisfaction.We highlighted that our design optimization approach resulted in a much smaller reservoir size with respect to the previous estimations (i.e., 315,000 vs. 1,100,000 m 3 ), with a minor loss of reliability.
The design procedure was then formalized by means of one-hundred Monte Carlo scenarios, in which we considered three key sources of uncertainty, involving physical (wind velocity), anthropogenic (load demand), and internal (wind power-to-electricity conversion) processes.As expected, the incorporation of all these sources of uncertainty within the design, which was the key novelty of this research, was reflected in all quantities of interest, namely, the reservoir size, the reliability of the system, and the net profit.Nevertheless, the initially proposed reservoir size still remained far away from the 90% range of variability.The quantification and interpretation of uncertainty was facilitated through typical statistical metrics (e.g., empirically derived quantiles), as well as more advanced probabilistic-stochastic tools.In this context, we applied a copula-based approach to evaluate the trade-offs between two of the most important aspects of the optimization problem, namely, the active depth of the storage component (design variable) and the net annual profit (performance metric).Our analysis indicated that this could also serve as a decision support tool for all the associated groups of interest (engineers, investors, and stakeholders).
Furthermore, the selection of a seawater PHS system introduced several challenges, which we briefly highlighted and addressed through effective technological means.To date, this option has only been practically tested in very limited cases; yet, it has the potential for wide applicability in islanding areas that suffer from water scarcity (which was also the case for Sifnos).Our analysis indicated that HRESs coupled with seawater PHS systems can provide techno-economically feasible solutions to the matter of energy independence in isolated and islandic regions, thus contributing to the sustainable development of the associated local communities.
Certainly, the overall design of HRESs supported by PHS systems also presented several other challenges to be addressed involving the project siting and synergies between the systems' components, thus making the underlying multi-objective optimization problem even more complex [95,96].
In particular, the generic layout of the system should also be investigated from the environmental perspective, since ecologically sensitive areas require careful planning and thorough assessing of the ecosystem impacts.Thus, an HRES feasibility study must be accompanied by a formalized environmental impact assessment (EIA).An EIA includes measures to minimize or offset the impacts of the surrounding and broader ecosystem, deriving from the project's implementation and operation [97].Mitigation measures refer to scaling down or even relocating the project to reduce the impacts at the source, whereas offset measures compensate for the negative impacts by providing solutions that counterbalance them [98].
On the other hand, the synergy between water and energy components across an HRES can be even more effective in order to improve its overall performance.A prime example is the installation of floating photovoltaic modules (FPVs) on a reservoir [99,100].Such a configuration can favor both components, since the water cooling effects of evaporation and wind ventilation increase the PV energy yield, while the PV modules' shading effects concurrently limit the losses due to evaporation [101].Such a layout can also be investigated in the Sifnos case.
The incorporation of environmental issues within the design optimization procedure and the interactions between the reservoir and FPVs introduced additional sources of uncertainty to be accounted for in a future study.There is also an ample number of uncertainty issues to be included in a forthcoming analysis with respect to external drivers (e.g., solar radiation, energy market, and social reactions), as well as power delivery issues (e.g., electricity transmission losses, voltage stability, and breakeven grid extension distance) [102].Nevertheless, this study proposed all the fundamental elements, in terms of the theoretical background and computational tools, to address the design of HRESs under uncertainty.

Figure 1 .
Figure 1.(a) The island of Sifnos and the HRES proposed siting (red circle) and (b) the proposed HRES layout (source: Google Earth, after processing).

Figure 2 .
Figure 2. Configuration of simulation model and associated water energy fluxes.

Figure 1 .
Figure 1.(a) The island of Sifnos and the HRES proposed siting (red circle) and (b) the proposed HRES layout (source: Google Earth, after processing).

Figure 1 .
Figure 1.(a) The island of Sifnos and the HRES proposed siting (red circle) and (b) the proposed HRES layout (source: Google Earth, after processing).

Figure 2 .
Figure 2. Configuration of simulation model and associated water energy fluxes.Figure 2. Configuration of simulation model and associated water energy fluxes.

Figure 2 .
Figure 2. Configuration of simulation model and associated water energy fluxes.Figure 2. Configuration of simulation model and associated water energy fluxes.

Figure 6 .
Figure 6.Fitting of the normal (left) and log-normal distributions (right) to the mean annual profit and reservoir active depth, respectively, where (a) displays the theoretical density functions and (b) displays the cumulative density functions.

Figure 6 .
Figure 6.Fitting of the normal (left) and log-normal distributions (right) to the mean annual profit and reservoir active depth, respectively, where (a) displays the theoretical density functions and (b) displays the cumulative density functions.

Figure 6 .
Figure 6.Fitting of the normal (left) and log-normal distributions (right) to the mean annual profit and reservoir active depth, respectively, where (a) displays the theoretical density functions and (b) displays the cumulative density functions.

Figure 7 .
Figure 7. Fitting of Gaussian copula to mean annual profit and reservoir active depth.

Figure 7 .
Figure 7. Fitting of Gaussian copula to mean annual profit and reservoir active depth.

Table 1 .
Wind turbine and solar PV properties.

Table 2 .
Unit economic data used in the simulation.

Table 4 .
Statistical characteristics of the optimized scenarios (one hundred datasets).

Table 4 .
Statistical characteristics of the optimized scenarios (one hundred datasets).