Mixing renewable energy with pumped hydropower storage: Design optimization under uncertainty and other challenges

1. Introduction

Hybrid renewable energy systems (HRES), which were first introduced in the 1970s, have become an essential part of global energy production, addressing limitations in terms of fuel flexibility, efficiency, reliability, emissions, and economics [1]. They generate electricity from two or more energy sources, while sharing a single connection point, allowing them to switch between energy sources when one is insufficient, thus reducing the inherited unpredictability of renewable energy sources. HRES can also capitalize on existing energy infrastructure and add components to help reduce costs, environmental impacts, and system disruptions [2]. Power can be generated based on the demand at any particular site depending on the availability of resources, thus also significantly reducing grid dependence.

The optimal design of an HRES 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 [3]. Interestingly, the absence of an 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 [4].

For this reason, the rapid development of hybrid renewable energy systems 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% [5,6,7,8]. Energy storage techniques are distinguished into four basic categories, according to their applications [9]. Low and medium-power applications (i.e., batteries, hydrogen fuel cells, superconductors) concern isolated areas and are used in emergency terminals or individual electrical systems. On the other hand, network connection with peak leveling and power-quality control applications (i.e., hydraulic systems, accumulators and flow batteries, compressed air) are used in power systems that provide energy for greater regions.

Pumped hydro storage (PHS) is the predominant and most reliable energy storage technology, accounting for more than 95% of global cumulative energy storage capacity [10]. Among its various benefits, Rehman et al. [11] highlighted several issues that PHS systems can address, which emerge from the large integration of wind power into the electricity network, such as: (1) handling changes in network impedances due to wind farm connection to the grid and its effects on the remote control-signals, (2) handling of harmonics created by the addition of wind on the grid and (3) stability problems that may occur due to dynamics behaviors 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 [12]. 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 been suggesting their economic feasibility, especially in remote areas and islandic systems [13,14,15,16].

HRES combined with PHS (which is in fact a water-energy system) are mainly applicable in non-interconnected power systems, such as the ones found in islands. These are defined by geographical isolation, which constitutes their underwater interconnection to the main grid economically not viable, low power demand, and high electricity production cost due to the importing of fossil fuels [17]. In particular, Greece could greatly benefit from the aforementioned systems, since it consists of 29 non-interconnected islandic complexes [18] and its renewable energy potential is high due to the favorable weather conditions. Moreover, the islands’ population increase during the summer leads to increased peak power demand, which further intensifies the need for energy independence. Today, only two HRES have been installed in non-interconnected Greek islands (i.e., Ikaria, Tilos), as pilot applications, while there are also few other systems under investigation, mainly in research context (e.g.,[19,20,21,22,23]). However, very few of these studies attempt to handle the major issue of uncertainty within the planning and design of such systems (e.g., [24]).

In general, PHS systems are configured as either open- or closed-loop. 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 PHS is part of a HRES, the storage capacity of the reservoir system is dictated by the power accounting between the highly fluctuating, 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 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, a HRES combining power production from renewables and energy storage through a PHS utilizing seawater was proposed for the island of Sifnos [17]. Its preliminary design followed the aforementioned conventional approach, resulting in a quite large reservoir capacity, if compared to the project scale. This case will serve as a proof of concept to further expand a recently proposed stochastic simulation-optimization framework to the design of HRES [25].

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 significantly 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. A parallel objective, originating from the use of the sea as lower reservoir, is to reveal the multiple technical challenges and resulting uncertainties, due to the use of seawater across multiple components of the PHS (corrosion of piping system, pumps and turbines, and risk of groundwater contamination in case of leakage).

2. 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² and a permanent population of 2,755 residents, as per the 2021 census conducted by the Hellenic Statistical Authority (ELSTAT). 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 a 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 (HEDNO), 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 [17] investigated the development of a 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 utilizes seawater. The upper reservoir is sited in a plateau, at an elevation of +320 m. This configuration is ideal, as there are no physical obstacles hindering the power production by the wind turbines and the PV’s. Moreover, the favoring topography (i.e., steep yet constant terrain slope) facilitates 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 [26], we highlight that the siting location is secluded, with the nearest settlement being 4 km away.

Key outcome of the preliminary design analysis was the sizing of the PHS, resulting to a quite large storage capacity of 1,100,000 m³. This capacity ensures 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 attempt to revise some design quantities, and particularly the reservoir size, following a more integrated simulation-optimization procedure (section 3). Furthermore, we embed the key sources of uncertainty within the design procedure, which is formalized in stochastic terms (Section 4).

3. HRES simulation and optimization

3.1. Configuration and key assumptions of 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 follows the one of [17], which has been further specified, after preliminary analysis. Its components are illustrated in the schematic sketch of Figure 2, while its key characteristics are summarized in Table 1. All essential design properties are specified, except for the reservoir active depth and the number of PV modules, which are considered subject to optimization.

The wind turbines’ tower height determines the wind speed values at the hub, which are logarithmically distributed with respect to 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 is expected that the interaction between large and small wind turbines (e.g., due to turbulence effects) will reduce the wind speed at the hub of the latter. The adjusted speed is calculated through the formula [27]:

$${\mathrm{V}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}}={\mathrm{V}}_0\left(1-\frac{2\mathrm{a}}{{\left(1+\frac{2\mathrm{k}\mathrm{L}}{{\mathrm{D}}_{\mathrm{L}}}\right)}^2}\right)$$

(1)

where V₀ 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 large turbine’s blade diameter. The a and k values are set equal to 0.10 and 0.038, respectively, as suggested in [27]. In our case, we consider the use of four wind turbines, i.e., two large of 2.3 MW and two small of 0.9 MW nominal power, respectively, distanced at 400 m. It is important to remark that greater distance might seem ideal, since the wind speed reduction is less, but does not necessarily result in optimal layouts [28].

The wind power production for a given wind speed value is calculated by the analytical formula, introduced by Sakki et al. [29]:

$$P_{wind}=P_{\mathit{min}}\begin{array}{c}+\left(1-{\left(1-{\left(\frac{V_{wind}-V_{min}}{V_{max}-V_{min}}\right)}^a\right)}^b\right){(P}_{max}-P_{min})\\ \end{array}$$

(2)

where P_min is the minimum power produced at 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 shape parameters. The two parameters are calibrated against empirically derived power curve data provided by the manufacturer. For the wind turbine types used in our simulation, we get a=2.25 and b=20. We remark that 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 ([30,31,32]).

Table 1. Wind turbine and solar PV properties.

	Wind turbines
Model	Enercon E-44	Enercon E-70 E4
Rated power (kW)	900.0	2,300.0
Minimum power (kW)	4.0	2.0
Cut-in wind speed (m/s)	3.0	2.5
Rated wind speed (m/s)	16.5	15.0
Cut-out wind speed (m/s)	34.0	34.0
Tower height (m)	55.0	113.0
Rotor diameter (m)	44.0	71.0
	Solar panels
Surface area (m2)	1.94
Nominal power (W)	340.0
Efficiency (%)	17.5

Table 1. Wind turbine and solar PV properties.

	Wind turbines
Model	Enercon E-44	Enercon E-70 E4
Rated power (kW)	900.0	2,300.0
Minimum power (kW)	4.0	2.0
Cut-in wind speed (m/s)	3.0	2.5
Rated wind speed (m/s)	16.5	15.0
Cut-out wind speed (m/s)	34.0	34.0
Tower height (m)	55.0	113.0
Rotor diameter (m)	44.0	71.0
	Solar panels
Surface area (m2)	1.94
Nominal power (W)	340.0
Efficiency (%)	17.5

The following assumptions are 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 at 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 is 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 during winter in the night;
• The turbine’s power capacity is also 6.0 MW, which is slightly larger than the maximum hourly load (5.4 MW) in order to account for uncertainties, as discussed later;
• The total efficiency values of turbines and pumps are considered constant and equal to 0.85 and 0.80, respectively;
• The penstock’s length and diameter are 910 m and 1.0 m, respectively, as specified in our preliminary design analysis.

3.2. Breakdown of simulation model

The simulation model represents the system’s operation in hourly time intervals, model and is explained through a flow chart in Figure 3. At each time step, the power accounting is performed, by contrasting the power production by the two renewables (wind, solar) to the actual demand. If there are energy surpluses (P_Net > 0), the PHS is set to its charging phase, thus pumping water from the sea to the upper reservoir, provided that there is sufficient storage capacity (S_total < S_max). Similarly, if there are energy deficits (P_Net < 0), the discharge phase starts and the water is released downstream through the turbine, thus generating electrical energy, provided that there is available water stored in the upper reservoir (S_total > V_dead). This process is repeated until the final step equals the simulation length, which in our case is 20 years. The entire simulation procedure is formulated in the Rstudio programming environment.

3.3. Setup of optimization problem

The 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 are often conflicting ([33,34]).

As already mentioned in Section 3.1, in our specific case we consider as design variables of the system the active depth of the reservoir and the number of PV modules. The optimization of these variables is formulated as a multi-objective problem, by accounting for different performance criteria that are expressed in financial terms and weighted properly, in order to ensure a well-balanced operation of the hybrid renewable energy system.

Specifically, the overall performance measure seeks for contrasting the revenues from electricity production, with the costs of construction, purchase, installation and maintenance (all expressed in terms of Equivalent Annual Cost, for the full depreciation time of the project, i.e., 20 years). The project implementation costs that were accounted for in the economic evaluation are associated with:

• the civil engineering works (excavations, roadworks, etc.);
• the purchase, installation and maintenance of the electromechanical equipment (wind turbines, PVs, pumps, turbines) and the conveyance system (GRP pipes);
• specific works associated with the reservoir waterproofing;

We remark that the aforementioned cost values in this study are following the typical pricing of respective projects in Greece.

A crucial requirement of electricity production systems is their reliability, which is a probabilistic concept, defined as the frequency of load demand satisfaction. In our approach, we utilize the complementary definition of reliability, i.e., failure probability, which is quantified in economic terms and embedded in the optimization by means of an empirical penalty function. In this respect, the optimization aims at ensuring a financially attractive and technically sound investment.

For the optimization, we utilize the evolutionary annealing-simplex algorithm [35]. EAS is a probabilistic heuristic global optimization technique, combining the robustness of simulated annealing in rough response surface with the efficiency of hill-climbing methods in convex areas.

3.4. 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³, and a PV power capacity of 1.09 MW. Furthermore, the key metrics of the optimized benchmark scenario are presented in Table 2. We remark that the proposed solution ensures 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 underline that the small capacity factor of the hydropower station (actual vs potential energy production) does not indicate a low performance. In contrast, it reveals 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 are concerned, they are in line with the climatic regime of the study area.

Interestingly, the optimized reservoir size equals less than a third of the one suggested in [17]. However, the two solutions are 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.

4. 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 HRES, deriving from multiple drivers. Sakki et al. [25] discriminated 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 utilizes the aforementioned stochastic optimization framework and applies it to the case of Sifnos, by embedding three key sources of uncertainty (i.e., wind process, energy demand, 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.

4.1. Wind process uncertainty

Hydrometeorological processes are considered one of the main exogenous uncertainties of a HRES due to the intermittent nature of renewables. Specifically, these concern wind processes (i.e., wind velocity, 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 is one of the key components of our uncertainty analysis, the mainstream approach to represent the variability of wind velocity is 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 in 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. [36] 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 is referred to as long-term persistence and is associated with the changing climate [37,38,39].

4.2. 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 [40]. Uncertainty in electricity demand, both for simulation and forecasting purposes, is a topic of high interest in literature. For instance, a study conducted by Cabeza et al. [41] showed that past energy projections for the energy demand of OECD countries were either consistently overestimated or underestimated. This phenomenon intensifies even more in the islands due to the intensive seasonal power demand variations. Thus, it is essential to accurately represent energy demand in probabilistic terms while designing a HRES. Warren et al. [42] also highlighted three main constituents of uncertainty in energy demand; (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. [43] presented various models for short, medium, and long-term energy demand modeling under uncertainties and metrics to evaluate their accuracy effectively.

4.3. Wind-to-power conversion uncertainty

As the penetration of wind power in energy systems increases, several concerns about the uncertainty in wind power generation are raised. Uncertainty in wind power system operations can be categorized between discrete and continuous disturbances [44]. The discrete disturbances to equipment failure, such as generators and transmission lines, were accurately demonstrated in an analysis performed on wind turbine generators by Rezamand et al. [45]. The results indicated 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 [46]. Given that wind turbines are commercial products, the WTPCs are provided by the manufacturers as a graph or as a set of points with a typical wind speed discretization of 0.5 m/s. Actually, these are theoretical relationships exhibiting a wind turbine's standard and experimental behavior [47]. 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 [48]. The deviation of wind turbine on-site behavior from the theoretical power curve is thoroughly analyzed in the work of Antoniou et al. [49], whilst several real-world examples are also presented in [50]. For this reason, an ample of deterministic and probabilistic models were developed to produce WTPCs that resemble actual operating conditions. Recently, the focus has been 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 [51]. A novel probabilistic WTPC model worth mentioning is the one developed by Yan et al. [52], which considers various model inputs (pitch angle, wind direction) based on three non-parametric algorithms.

A last major issue of uncertainty involves the operation of wind turbines in the high wind speed 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. However, due to the stochastic nature of wind, and the uncertainty of the theoretical WTPCs, it is impossible to model in deterministic terms when the wind turbines will shut down, and their downtime [53]. Thus, to prevent frequent shutdowns and restarts, the soft cut-out strategy is implemented by extending the maximum admissible wind speed up to 30-32 m/s, without an abrupt shutdown, through controlling pitch and generator in order to gradually decrease the energy production [54]. 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. [55] produced a soft-cut out strategy worst-case scenario algorithm that does not significantly increase fatigue loads. Astolfi et al. [56] 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. [54] 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 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, is associated with monitoring errors. The output power observations are 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 [57]. Consequently, these data need to be cleared by appropriate classification algorithms in order to eliminate their adverse effects and improve wind power predictability [58].

5. HRES simulation and optimization under uncertainty

5.1. Incorporating uncertainty in the simulation

Following the uncertainty-aware framework by Sakki et al. [25], we run the design optimization problem in a Monte Carlo context. This involves the formulation of a large number of scenarios (in the specific case, one hundred), which allows 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 are embedded in the aforementioned Monte Carlo simulation-optimization approach, distinguishing 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 are integrated in the simulation by means of synthetic time series that are generated through the anySim package [39], which provides a suite of stochastic models for the simulation of both stationary and cyclostationary processes (in univariate or multivariate mode) that may follow a wide spectrum of distributions. This allows for generating synthetic time series with the desired marginal and stochastic properties, as reflected in the associated historical data [59,60,61]. In our work, we utilize this package, which runs 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 HRES. Their probabilistic and dependence regime reproduces the one of the corresponding historical data, the length of which (ten years) is rather small for a proper representation of their actual variability, and thus uncertainty.

Regarding the uncertainty in 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 do not accurately represent their on-site behavior. In this vein, we follow a Monte Carlo simulation approach, considering that the shape parameters a and b of the empirical Equation (2) are random variables. Specifically, for both turbines the output power for a given wind velocity value is sampled from a normal distribution with a ~ N(2.25, 0.0225) and b ~ N(20, 0.016). This distribution and the associated parameters (i.e., mean and standard deviation values) were appropriately selected, after some trials, so that the average deviation from the theoretical curve is no greater than 15% (Figure 4). This deviation is within reasonable margins, since Lira et al. [62] stated that power curve uncertainty is approximately 10%.

Furthermore, we also implement the soft cut-out strategy, allowing the wind turbine to operate in the high-wind speed region (Figure 4). Therefore, we consider a linear reduction of the turbine’s power for wind speed values between 25 and 30 m/s. The slope, z, of this linear equation is also represented as a normal process, with z ~ N(150, 225) for the large turbines and z ~ N(45, 100) for the small ones.

5.2. Results of Monte Carlo scenarios

Table 3 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 remind that the two design variables to optimize were the reservoir active depth and the total power capacity of 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, exhibits significant variability, which is propagated in the key performance metrics of the system (net profit and reliability). In contrast, the variability of PV capacity is negligible. In this context, the range of uncertainty with respect to renewable energy production mainly reflects the stochasticity of the wind process and its conversion through the probabilistic power curves of the two turbines. We remark that while the installed wind power capacity remains constant (i.e., 6.4 MW), the capacity factor of the two machines varies (although not significantly), as result of this stochasticity. On the other hand, since for all scenarios we apply the same set of solar radiation data, the capacity factor of PVs remains constant, and equals 0.21.

5.3. Insight into the trade-off between reservoir size and overall system profit

The analysis of the model outcomes under uncertainty reveals the existence of a trade-off between the key design quantity, i.e., the active depth of the reservoir, which dictates its size, and the overall profit of the water-energy system. Since this trade-off embeds all kinds of uncertainties that are included in the simulation-optimization procedure, its quantification and interpretation presuppose the use of advanced stochastic and statistical tools. This objective is implemented via a two-step approach.

First, we investigate the marginal behavior of the two random variables, by assigning suitable statistical distributions, specifically the normal and the log-normal one, to the reservoir active depth and the mean annual profit, respectively (Figure 5). In this context, we utilize the fitdistrplus R package [63], which provides 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 apply the moment-matching estimation method, which involves finding the values of the model parameters that make the data's sample moments equal to the model's corresponding population moments.

After selecting the marginal distributions for the two variables, we investigate their dependencies, by fitting a Gaussian copula, as shown in Figure 6. Notably, the two variables exhibit a quite strong negative correlation, as indicated by the associated Pearson coefficient (R), which equals -0.75. The copula approach enables 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 modelling of dependent variables of any type ([64,65]), as well as for predictive uncertainty assessment studies.

In our case, the copula-based approach provides further 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 perspective of investors and stakeholders.

6. The challenge of seawater

One of the major obstacles for the development of HRES in small, non-interconnected islands, which take advantage of PHS, is the scarcity of water resources, which is a long-lasting issue, especially during summer. A potential solution could be the deployment of PHS systems that utilize the sea as the lower reservoir, as proposed in the case of Sifnos island. The Okinawa Yanbaru Seawater Pumped Storage Power Station [66] was the first experimental pumped-storage facility worldwide that used seawater for energy storage. Its construction lasted from 1987 to 1999, yet it was dismantled in 2016. This is because Okinawa’s projections overestimated the load demand growth, constituting the plant’s operation no longer profitable.

Over the years, there have been several studies suggesting the use of PHS systems that utilize the sea as the lower reservoir in islands, including Hawaii, Ireland and Azores [67,68,69]. In a more recent study, Ali and Jang [70] examined the case of Deokjeok-do, South Korea and found that a HRES coupled with a seawater PHS, apart from providing energy independence, also leads to a cheaper initial investment cost than using batteries as energy storage. However, the inclusion of a seawater pumped storage in the HRES introduces various technical challenges, mainly related to corrosion. This section aims to address the most crucial ones, ensuring the proposed system’s technical feasibility.

6.1. Conveyance system

Corrosion prevention can vary depending on whether the components are subject to high or low flow. This protection categorization was implemented in the Okinawa project [71]. 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) is chosen for smaller pipe diameters [72]. 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 [73]. 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 [72].

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 [74]. 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 [75]. However, it was found that such practices not only threatened the surrounding ecosystem, but also deposited byproducts (i.e., manganese) [76]. While less harmful methods that do not produce biocidal intermediates have been investigated, i.e., electrochemical degradation ([77]), natural biocide-based and non-stick coatings ([78,79]), further research is required in order to provide solutions that ensure both the surrounding ecosystems’ sustainability and the projects’ techno-economic feasibility.

6.2. Electromechanical equipment

Electromechanical equipment refers to the turbines and pumps used in the PHS system. In HRES, the selection of pumps is a particularly important procedure, as there is a wide variety available in 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 [80]. 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 head to fill the reservoir. This issue was addressed by Manfrida and Secchi [80], 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 [81] distinguished the several types of corrosion from seawater that can affect stainless steel (SS) electromechanical equipment into the following:

• 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 tend to occur in crevices. In a seawater pumps, crevices can be found where seals and impellers are fastened to the shaft, and flange faces are cast in for pipe-work 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, Molybdenum) have a high pitting resistance equivalent number (PREN), which indicates its corrosion resistance. According to the Norwegian Standards (NORSOK) a PREN number of 40 provides sufficient corrosion resistance to seawater [82].

6.3. 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 by wind transport). Specifically, groundwater salinization is a crucial environmental issue, especially in areas with karstified underlying geological formations (this is 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 (HDPE), which are resistant to chemicals and UV radiation, are a potential solution, combined with a drainage system that detects possible leakages [66]. On the other hand, the construction of an embankment around the reservoir can prevent seawater transport through the wind [72].

7. Conclusions

In the present study, we developed a simulation-optimization procedure for the design of HRES 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 fluxes. With regards 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 load demand satisfaction.

The design procedure was formalized by means of one hundred Monte Carlo scenarios, in which we considered three key sources of uncertainty, involving a physical process (wind velocity), an anthropogenic process (load demand) and an internal process (wind power-to-electricity conversion). 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 can also serve as a decision support tool for all associated groups of interest (engineers, investors, stakeholders).

Furthermore, the selection of a seawater PHS introduces several challenges, which we briefly highlighted and addressed through effective technological means. To this day, this option has only been practically tested in very limited cases, yet has the potential for wide applicability in islanding areas that suffer from water scarcity (which is also the case for Sifnos). Our analysis indicated that HRES coupled with seawater PHS 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.

Certainty, the overall design of HRES supported by PHS also comprise several other challenges to address, involving the project siting and the synergies between the systems’ components, thus constituting the underlying multiobjective optimization problem even more complex ([83,84]).

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, a HRES feasibility study must be accompanied by a formalized environmental impact assessment (EIA). EIA includes measures to minimize or offset the impacts of the surrounding and broader ecosystem, deriving from the project’s implementation and operation [85]. 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 [86].

On the other hand, the synergy between water and energy components across a HRES can be even more effective, in order to improve its overall performance. A prime example is the installation of floating photovoltaic modules (FPV) on the reservoir [87,88]. Such a configuration can favor both components, since the water cooling effects of evaporation and wind ventilation increase PV energy yield, while the PV modules’ shading effects concurrently limits the losses due to evaporation [89]. Such a layout could also be investigated for the Sifnos case.

The incorporation of environmental issues within the design optimization procedure, and the interactions between the reservoir and the FPVs introduce additional sources of uncertainty, to be accounted for in a future study. There is also an ample of uncertain variables to be included in a forthcoming analysis with respect to external drivers (e.g., solar radiation, energy market) as well as power delivery (e.g., electricity transmission losses, voltage stability and breakeven grid extension distance) [90]. Nevertheless, this study has set all fundamental elements, in terms of theoretical background and computational tools, to address the design of HRES under uncertainty.