Transient Simulation of Underground Pumped Storage Hydropower Plants Operating in Pumping Mode

: The increasing penetration of variable renewable energies (VRE) in the European electricity mix requires ﬂexible energy storage systems (ESS), such as pumped storage hydropower (PSH). Disused mining voids from deep closed mines may be used as subsurface reservoirs of underground pumped-storage hydropower (UPSH) plants. Unlike conventional PSH plants, the air pressure in UPSH plants is variable and it di ﬀ ers from the atmospheric conditions. In this paper, the hydraulic transient process of an UPSH plant operating in pumping mode was investigated and a preliminary thermodynamic analysis of the closed surge tank was carried out. Analytical and CFD three-dimensional numerical simulations based on the volume of ﬂuid (VOF) model with two-phase ﬂow have been performed for analyzing the transient process. In the transient simulation, air and water are considered as ideal gas and compressible liquid, respectively. Di ﬀ erent guide vanes closing schemes have been simulated. The obtained results show that the dimensioning of underground reservoir, surge tank, and air ducts is essential for ensuring the hydraulic performance and optimizing the operation of UPSH plants. The static pressure in the air duct, surge tank and lower reservoir reaches − 1.6, 112.8 and − 4 kPa, respectively, while a heat ﬂux of − 80 W was obtained through the surge tank walls. heat flux also increase, varying the thermodynamic behavior. This study demonstrates the paramount importance of the dimensioning of the underground infrastructure of UPSH plants. The obtained results show that the thermodynamic behavior during transient processes depends on the surge tank size, rock mass material, and period of the mass oscillations.


Introduction
In 2019, renewable energies generated 34.6% of Europe's electricity [1]. Renewable generation (excl. hydro) increased by 65 TWh in 2019, exceeding the 2010-2018 average of 50 TWh per year. This high level of renewable energy sources (RES) production was mainly attributable to an increase in wind generation. Since 2015, 84% of the growth in total renewables generation has come from wind, 18% from solar, and 10% from biomass [1]. The increasing of some ways of intermittent renewable energy sources (RES), such as wind and solar PV, is making energy storage systems (ESS) increasingly important for balancing electricity supply and demand [2,3]. Pumped-storage hydropower (PSH) plants, as the most efficient and mature large energy storage facilities, will play an essential role in the power grid [4,5]. However, topographical requirements (elevation difference between the two reservoirs as well as large water volumes) and environmental concerns reduce the viable sites for the construction of PSH plants [6,7]. An alternative to solving this inconvenience might be the use of disused underground structures, such as deep closed mines. Underground pumped-storage hydropower (UPSH) uses an upper reservoir that provides water storage capacity at ground level, and a lower reservoir at a depth underground to give a suitable pressure. The powerhouse itself, accommodating Francis pump-turbine and motor-generator units, surge tank, and air ducts are also underground. The main advantage of UPSH plants is that suitable upload topography is not required and the environmental impact could be minimized. The idea of locating the lower water reservoir of a PSH scheme underground was proposed in 1901 by Reginald A. Fessenden [8]. Some studies of UPSH plants have been developed in recent years [2,6,[9][10][11][12][13][14]. However, there are no bibliographic evidence of UPSH plants under operation worldwide.
The operation of the UPSH plants is complex due to the presence of water and air interacting in the operation phase. The air pressure influences the energy production and efficiency of UPSH plants. Precisely, the authors have previously studied the effect of the air pressure on the global efficiency of the UPSH plants [15]. The round trip efficiency could be reduced down to 5% when the reservoir pressure reaches −100 kPa. Moreover, a study for determining the effect of the air pressure on the energy production of UPSH plants operating in turbine mode was also developed by Menendez et al. [16]. The reduction in the energy generation could reached up to 62 MWh cycle −1 (12.5% of the available energy) when the cross section of air ducts was too small. The air pressure inside the lower reservoir and air ducts depends on the water flow rate and the cross section of air ducts. The net head also depends on the air pressure inside the lower reservoir in UPSH plants, unlike conventional PSH plants.
Pummer et al. [17] analyzed the reflection phenomena in underground pumped storage reservoirs. The results shown excessive wave heights through wave reflections. The optimization of guide vanes closing schemes is an important strategy for reducing water hammer pressure and pulsating pressures [18]. The existence of this kind of pressure is one of the main reasons for unit vibration as well as fatigue damage of Francis pump-turbines [19,20]. Therefore, minimizing water hammer pressure during hydraulic transient processes is of great importance in ensuring the safety and stability of PSH plants. Vakil [21] examined different guide vanes closing laws to investigate their effects on pressure rise and speed rise. Zeng et al. [22] theoretically analyzed the effects of the guide vanes closing laws schemes on water hammer and pulsating pressures that are based on the transient characteristics of pump turbines in the S-shaped region. The hydraulic and thermodynamic behavior in the surge tank are primordial in transient processes. Thermodynamic behavior in closed surge tanks is usually studied with the polytropic equation, where a polytropic exponent value (Υ) is considered. [23]. The solution to such a process was the isothermal (Υ = 1) or isentropic assumption (Υ = 1.4) [24], although an average value Υ = 1.2 could be used in preliminary design calculations [25][26][27][28]. For fast transients, field observations and experiments show that the polytropic exponent value is approximately 1.4 and the thermodynamic behavior is close to adiabatic [29,30]. However, the behavior of the closed surge tanks for slow transients show that the heat transfer has a significant effect on the thermodynamic of the system and the polytropic equation does not properly represent the heat transfer. The water level fluctuation and air pressure depends on the dimensioning of the surge tank, the water flow rate, and the closing time of the guide vanes. Vereide et al. [31] carried out a hydraulic scale model of mass oscillations in closed surge tanks constructed as an underground rock cavern. The results revealed adiabatic behavior of the closed surge tank. Zhou et al. [32] studied the power-off transient characteristics of a PSH station in pump mode using a three-dimensional (3D) unsteady numerical method based on a single-phase and VOF coupled model.
In this work, the hydraulic transient process of an UPSH plant operating in pumping mode is analyzed. The air pressure inside the upper reservoir, surge tank, and air ducts is analyzed, depending on the dimensioning of the UPSH scheme. Different guide vanes closing schemes of a pump-turbine unit, modelled with user defined functions (UDF), have been introduced in the simulations to investigate their effects on pressure rise in air duct, surge tank, and lower reservoir. In addition, the thermodynamic behavior of the closed surge tank excavated in a sandstone rock mass Energies 2020, 13, 1781 3 of 17 is also studied. Internal energy, air temperature, and heat flux through the walls of the surge tank have been analyzed. Analytical and 3D numerical models that are based on volume of fluid (VOF) model with two-phase flow have been developed to analyze the operation of an UPSH plant during transient processes. The main objective is to highlight the paramount importance of considering the dimensioning of the underground reservoir, surge tanks, and air ducts to ensure the hydraulic performance and optimize the operation of UPSH plants.

Problem Statement
The problem is formulated, as shown in Figure 1. An UPSH plant in a closed coal mine at 450 m depth is considered. Table 1 shows the main geometrical data from the UPSH plant and the pump-turbine specifications. The flow rate at full load in pumping mode is 40 m 3 s −1 . The upper reservoir provides a water flow rate that passes through the Francis turbine generating electricity. Conversely, when the water pass from the lower reservoir to the upper reservoir, the system consumes energy (pumping mode). The round trip energy efficiency of UPSH plants (the energy fed into the electrical grid divided by the energy taken from the electrical grid in pumping mode), is in the 70-80% range [33,34]. The amount of storable energy depends on the water mass and net head between upper and lower reservoir [35]. In UPSH plants, the air flow direction through the ventilation shafts switches, depending on the operation mode. The air flow rate is constant in the pumping (downward direction) and turbine modes (upward direction), and it varies during hydraulic transient processes. during transient processes. The main objective is to highlight the paramount importance of considering the dimensioning of the underground reservoir, surge tanks, and air ducts to ensure the hydraulic performance and optimize the operation of UPSH plants.

Problem statement
The problem is formulated, as shown in Figure 1. An UPSH plant in a closed coal mine at 450 m depth is considered. Table 1 shows the main geometrical data from the UPSH plant and the pump-turbine specifications. The flow rate at full load in pumping mode is 40 m 3 s -1 . The upper reservoir provides a water flow rate that passes through the Francis turbine generating electricity. Conversely, when the water pass from the lower reservoir to the upper reservoir, the system consumes energy (pumping mode). The round trip energy efficiency of UPSH plants (the energy fed into the electrical grid divided by the energy taken from the electrical grid in pumping mode), is in the 70-80% range [33,34]. The amount of storable energy depends on the water mass and net head between upper and lower reservoir [35]. In UPSH plants, the air flow direction through the ventilation shafts switches, depending on the operation mode. The air flow rate is constant in the pumping (downward direction) and turbine modes (upward direction), and it varies during hydraulic transient processes.

Mathematical Formulation
A one-dimensional (1D) mathematical formulation has been developed for simulating the hydraulic transient process of an UPSH plant when considering different dimensioning of the surge tank and closing times of guide vanes (1, 10, and 15 s). In the present study, the fourth order Runge-Kutta method is used within the MATLAB software to solve the pressure drop oscillation. Therefore, a system of equations is applied iteratively, with a time step of 0.1 s. The energy equation between the underground reservoir (R) and the surge tank (S) in non-stationary regime has been applied in Equation (1) to estimate the variation of the water velocity inside the tunnel (see Figure 1). The direction of water flow rate inside the surge tank and tunnel varies during the transient process. Therefore, the absolute value of the water velocity and water level are considered in Equation (2). In Equation (3), the mass conservation is applied in the connection between the surge tank and the tunnel. The loss coefficient (ξ e ) in the expansion inside the surge tank is assumed to be ξ e = 0.5 when v > 0, and ξ e = 1 when the v < 0.

Mesh Description
A three-dimensional CFD numerical model is established to compare the results with the 1D analytical model. This numerical model is developed meshing the geometry into different small cells in which the flow equations must be solved after discretization. The meshes of the underground reservoir, tunnels, closed surge tank, and air duct were generated using the commercial software GAMBIT. A 3D structured mesh has been developed for the entire domain, employing tetrahedral and hexahedral cells because of their suitable adaptation to the geometry. Finally, the mesh was refined in the areas of interest, air duct, and closed surge tank. Throughout the calculations, the mesh for the underground reservoir, surge tank, and air duct was accurate with the maximum skewness less than 0.75. The mesh used for the calculations reaches up to about 2.74 × 10 6 cells and the initial number of nodes is 1.74 × 10 6 . Figure 2 shows the geometry of the considered UPSH scheme and the mesh details.  Figure 3 shows the guide vanes linear closing schemes. Figure 3a indicates the variation of the velocity inlet (-1.33 m s -1 in pumping mode at full load), and Figure 3b shows the guide vanes opening for the closing schemes considered. Water and air have been considered as compressible liquid and ideal gas, respectively. Table 3 illustrates the properties considered for water and air [41]. In addition, a Bulk Modulus (BM) for the water of 2.2 E+09 was also considered.

Mesh Grid Sensitivity
The mesh grid sensitivity analysis was carried out in order to evaluate the convergence of numerical solutions [36]. The performances of the CFD numerical models are heavily affected by the choice of the mesh size [37,38]. To this purpose, four simulations were carried out from t = 0 s to t = 5 s (operating in pumping mode) based on identical operating conditions and material properties. Table 2 compares the error on static pressure in the ventilation air shaft with four different mesh sizes. In addition, the computational time is also presented for each scenario. When the total grid cells is above 2.74 × 10 6 , the results are affected very little by the computational grid. A cluster with six cores and 64 GB of RAM was employed to perform the numerical simulations.

Solution Setup and Boundary Conditions
The objective of the simulation is to simulate the hydraulic transient process of an UPSH plant to obtain the values of the air pressure for different closing times inside the lower reservoir, air duct, and closed surge tank during transient processes. Water level fluctuation inside the surge tank has Energies 2020, 13, 1781 6 of 17 also been simulated. Finally, the thermodynamic behavior inside the surge tank has been analyzed during the transient process. Simulations were carried out with the pump-turbine specifications that are indicated in Table 1, while considering a single air duct with 1 m in diameter. The simulations were performed with the commercial CFD software Ansys Fluent V16.0, solving the Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations with a two-phase (water and air) scheme. The coupled momentum and continuity equations have been solved [39].
where ρ is density, v is velocity, P is pressure, S is stress tensor, and F is the contribution from external forces. The momentum Equation (4) expresses Newton's second law of motion, and the continuity Equation (5) ensures the conservation of mass.
In this model, a second-order upwind discretization has been implemented for convective terms, with a first order scheme for temporal terms and central differencing for diffusion terms. The VOF approach has been selected for the two-phase modelling. Its formulation is designed for two or more immiscible fluids where the position of the interface between the fluids is of interest [40]. The variables and properties of each phase are volume-averaged values, as calculated using a volume fraction ratio of each phase for every element. The model includes continuity equations for each phase and the URANS momentum equation for the entire domain. Those equations allow for describing each fluid movement as well as tracking the water liquid-air interface (considering both water and air as compressible flows). In a two-phase system (water and air), if the volume fraction of the second of these is being tracked, the density in each cell is given by Equation (6). The evolution of scalar f (volume fraction) is governed by the simple advection Equation (7) [40].
For the closure of turbulence, a robust k-epsilon Re-Normalisation Group (RNG) model with standard logarithmic wall functions was employed. A set of boundary conditions was selected to solve the equations in the transient process in pumping mode. Velocity inlet boundary condition has been employed at the entrance of the model (Figure 2), where different UDF were applied, while considering 1, 10, and 15 s of closure time of the guide vanes. Figure 3 shows the guide vanes linear closing schemes. Figure 3a indicates the variation of the velocity inlet (−1.33 m s −1 in pumping mode at full load), and Figure 3b shows the guide vanes opening for the closing schemes considered. Water and air have been considered as compressible liquid and ideal gas, respectively. Table 3 illustrates the properties considered for water and air [41]. In addition, a Bulk Modulus (BM) for the water of 2.2 × 10 9 was also considered. Table 3. Water and air properties considered in the numerical simulations.

Properties
Water Air   Pressure outlet boundary condition was considered at the exit of the ventilation shaft, while assuming a stable pressure value that corresponds to the standard atmospheric pressure (101.3 kPa). A surface tension coefficient of 0.073 N m -1 between water and air and a specified operating density of 1.225 kg m 3 were also selected. Finally, the unsteady resolution of the discretized equations, with time steps in the order of 10E-03 s, provides the final solution of the model. More precisely, variable time stepping was used to automatically change the time step when an interface moves through finer meshes or if the interface velocity is high. The variable time step is based on the maximum Courant number near the VOF interface.
The transient form of the Pressure-Implicit with Splitting of Operators (PISO) algorithm has been used to resolve the coupling between pressure and velocity fields in the Navier-Stokes equations. The spatial and temporal derivatives of the fluid flow governing equations were calculated by means of a second-order discretization, and, for the pressure interpolation, the PREssure STaggering Option (PRESTO) scheme has been used instead of the standard one.
The energy equation was applied in the models to analyze the thermodynamic behavior inside the closed surge tank in the transient process. The heat flux through walls of the surge tank the temperature distribution, enthalpy, and internal energy were estimated in the simulations. The closed surge tank is excavated in a sandstone rock mass. A specific heat of 711 J Kg -1 K -1 and a thermal conductivity of 5 W m -1 K -1 were considered in the rock mass at an initial temperature of 300 K (at 450 m depth). In the simulations, heat transfer is assumed to be dominated by the combination of convection in the air (closed surge tank) and conduction through the rock mass.

Analytical model results
The mathematical formulation has been applied for solving the hydraulic transient processes from t=0 s (beginning) to t=100 s (end). Different guide vanes closing schemes have been considered (1, 10, and 15 s, which correspond to the thick black, solid black, and dashed grey lines in Figure 4, respectively). Figure 4 indicates the surge tank air pressure rise. Adiabatic and isothermal behavior Pressure outlet boundary condition was considered at the exit of the ventilation shaft, while assuming a stable pressure value that corresponds to the standard atmospheric pressure (101.3 kPa). A surface tension coefficient of 0.073 N m −1 between water and air and a specified operating density of 1.225 kg m 3 were also selected. Finally, the unsteady resolution of the discretized equations, with time steps in the order of 10 × 10 −3 s, provides the final solution of the model. More precisely, variable time stepping was used to automatically change the time step when an interface moves through finer meshes or if the interface velocity is high. The variable time step is based on the maximum Courant number near the VOF interface.
The transient form of the Pressure-Implicit with Splitting of Operators (PISO) algorithm has been used to resolve the coupling between pressure and velocity fields in the Navier-Stokes equations. The spatial and temporal derivatives of the fluid flow governing equations were calculated by means of a second-order discretization, and, for the pressure interpolation, the Pressure Staggering Option (PRESTO) scheme has been used instead of the standard one.
The energy equation was applied in the models to analyze the thermodynamic behavior inside the closed surge tank in the transient process. The heat flux through walls of the surge tank the temperature distribution, enthalpy, and internal energy were estimated in the simulations. The closed surge tank is excavated in a sandstone rock mass. A specific heat of 711 J Kg −1 K −1 and a thermal conductivity of 5 W m −1 K −1 were considered in the rock mass at an initial temperature of 300 K (at 450 m depth). In the simulations, heat transfer is assumed to be dominated by the combination of convection in the air (closed surge tank) and conduction through the rock mass.

Analytical Model Results
The mathematical formulation has been applied for solving the hydraulic transient processes from t = 0 s (beginning) to t = 100 s (end). Different guide vanes closing schemes have been considered (1, 10, and 15 s, which correspond to the thick black, solid black, and dashed grey lines in Figure 4, respectively). Figure 4 indicates the surge tank air pressure rise. Adiabatic and isothermal behavior have been considered in Figures 4a and 4b      The surge tank water level fluctuation and the surge tank inlet velocity between t=0 s and t=100 s are shown in Figure 5a and Figure 5b, respectively. At the beginning of the closing process (t=0 s), the water level is -173 m. The maximum variation of the water level is 6.5 m (-166.5 m) at the instance t=6 s. The maximum water level is obtained when considering a closure law with duration of 1 s. The maximum water level fluctuation for tc=10 s and tc=15 s is 5 m and 3.5 m, respectively.         When the diameter and the height of the surge tank are reduced, the initial volume of air is reduced and, therefore, the pressure and the level of water are also reduced, extending the duration of the transient process. The maximum and minimum values of the air pressure in the closed surge tank are shown in Figure 7a and Figure 7b, respectively.

Numerical Results
Numerical simulations have been carried out from t = 0 until t = 40 s. The transient processes start at the instance t = 5 s. From the instance t = 0 s to t = 5 s, the pump-turbine is operating in pumping mode with the guide vanes opened at 100%. The initial water flow rate is 40 m 3 s −1 , so the water velocity in the tunnel is 1.33 m s −1 . In the pumping operation mode, the water level and the static pressure in the closed surge tank decrease and the static pressure in the air duct is constant. At the beginning of the closing process (t = 5 s), guide vanes linear closing schemes with a total duration of 1, 10 and 15 s have been applied. Figure 8 shows the water level at t = 5 s, t = 11.5 s, t = 14.5 s, and t = 19 s for a closure time of 1 s. Numerical simulations have been carried out from t=0 until t=40 s. The transient processes start at the instance t= 5 s. From the instance t=0 s to t=5 s, the pump-turbine is operating in pumping mode with the guide vanes opened at 100%. The initial water flow rate is 40 m 3 s -1 , so the water velocity in the tunnel is 1.33 m s -1 . In the pumping operation mode, the water level and the static pressure in the closed surge tank decrease and the static pressure in the air duct is constant. At the beginning of the closing process (t=5 s), guide vanes linear closing schemes with a total duration of 1, 10 and 15 s have been applied. Figure 8 shows the water level at t=5 s, t=11.5 s, t=14.5 s, and t=19 s for a closure time of 1 s.  Figure 9a shows the evolution of the air velocity in the ventilation shaft. The maximum value of the velocity is reached from t=0 to t=5 s (pumping mode). Note that the negative values correspond to a situation of inlet flow rate (inhalation from the atmosphere). Following, the air velocity magnitude decreases from the beginning of the closing process, at the instance t=5 s. Complementarily, the maximum static pressure in the air duct is -1.66 kPa in relative terms ( Figure  9b).    Figure 9a shows the evolution of the air velocity in the ventilation shaft. The maximum value of the velocity is reached from t = 0 to t = 5 s (pumping mode). Note that the negative values correspond to a situation of inlet flow rate (inhalation from the atmosphere). Following, the air velocity magnitude decreases from the beginning of the closing process, at the instance t = 5 s. Complementarily, the maximum static pressure in the air duct is −1.66 kPa in relative terms (Figure 9b). Numerical simulations have been carried out from t=0 until t=40 s. The transient processes start at the instance t= 5 s. From the instance t=0 s to t=5 s, the pump-turbine is operating in pumping mode with the guide vanes opened at 100%. The initial water flow rate is 40 m 3 s -1 , so the water velocity in the tunnel is 1.33 m s -1 . In the pumping operation mode, the water level and the static pressure in the closed surge tank decrease and the static pressure in the air duct is constant. At the beginning of the closing process (t=5 s), guide vanes linear closing schemes with a total duration of 1, 10 and 15 s have been applied. Figure 8 shows the water level at t=5 s, t=11.5 s, t=14.5 s, and t=19 s for a closure time of 1 s.  Figure 9a shows the evolution of the air velocity in the ventilation shaft. The maximum value of the velocity is reached from t=0 to t=5 s (pumping mode). Note that the negative values correspond to a situation of inlet flow rate (inhalation from the atmosphere). Following, the air velocity magnitude decreases from the beginning of the closing process, at the instance t=5 s. Complementarily, the maximum static pressure in the air duct is -1.66 kPa in relative terms ( Figure  9b).    The static pressure in the longitudinal section of the model might be observed in Figure 11 at the beginning of the closing process ( Figure 11a) and at the instance t=16.5 s (Figure 11b) when considering a closing scheme of guide vanes with a duration of 10 s. Static pressure values in the air duct, lower reservoir, tunnel, and surge tank are shown. At the beginning of the guide vanes closing process, specifically at time t=5 s, the air manometric pressure in the surge tank reaches 40 kPa. The air pressure increases up to 90 kPa at time t=16.5 s. Obviously, the static pressure increases up to 320 kPa in the tunnel due to the water column. The static pressure in the longitudinal section of the model might be observed in Figure 11 at the beginning of the closing process ( Figure 11a) and at the instance t = 16.5 s (Figure 11b) when considering a closing scheme of guide vanes with a duration of 10 s. Static pressure values in the air duct, lower reservoir, tunnel, and surge tank are shown. At the beginning of the guide vanes closing process, specifically at time t = 5 s, the air manometric pressure in the surge tank reaches 40 kPa. The air pressure increases up to 90 kPa at time t = 16.5 s. Obviously, the static pressure increases up to 320 kPa in the tunnel due to the water column. A preliminary thermodynamic analysis has been carried out in the transient processes. Figure  12 depicts the variation of the internal energy (U) and the evolution of the air temperature for the three closing schemes considered. The internal energy and the air temperature increase when the pressure raises. The variation of the internal energy reaches 3,450 J kg -1 K -1 at the instance t=6 s for a closing scheme with a duration of 1 s. When the static pressure is reduced, the variation of the internal energy is also reduced down to -3,008 J kg -1 K -1 at t=15 s. The variation of internal energy is lower for closure times of 10 s and 15 s. The maximum value of air temperature is 333 K for a closure time of guide vanes of 1 s. This is obtained when the air pressure reaches 112.8 kPa at the instance t=11.5 s. The air temperature for the closure times of guide vanes of 10 s and 15 s is lower, reaching 324 K and 315 K, respectively. A preliminary thermodynamic analysis has been carried out in the transient processes. Figure 12 depicts the variation of the internal energy (U) and the evolution of the air temperature for the three closing schemes considered. A preliminary thermodynamic analysis has been carried out in the transient processes. Figure  12 depicts the variation of the internal energy (U) and the evolution of the air temperature for the three closing schemes considered. The internal energy and the air temperature increase when the pressure raises. The variation of the internal energy reaches 3,450 J kg -1 K -1 at the instance t=6 s for a closing scheme with a duration of 1 s. When the static pressure is reduced, the variation of the internal energy is also reduced down to -3,008 J kg -1 K -1 at t=15 s. The variation of internal energy is lower for closure times of 10 s and 15 s. The maximum value of air temperature is 333 K for a closure time of guide vanes of 1 s. This is obtained when the air pressure reaches 112.8 kPa at the instance t=11.5 s. The air temperature for the closure times of guide vanes of 10 s and 15 s is lower, reaching 324 K and 315 K, respectively. The internal energy and the air temperature increase when the pressure raises. The variation of the internal energy reaches 3450 J kg −1 K −1 at the instance t = 6 s for a closing scheme with a duration of 1 s. When the static pressure is reduced, the variation of the internal energy is also reduced down to −3008 J kg −1 K −1 at t = 15 s. The variation of internal energy is lower for closure times of 10 s and 15 s. The maximum value of air temperature is 333 K for a closure time of guide vanes of 1 s. This is obtained when the air pressure reaches 112.8 kPa at the instance t = 11.5 s. The air temperature for the closure times of guide vanes of 10 s and 15 s is lower, reaching 324 K and 315 K, respectively. Figure 13 shows the total temperature and the heat flux through the walls of the closed surge tank for both water and air phases. The water temperature remains constant throughout the transient process. However, the air temperature increases when the static pressure increases. A heat transfer of −80 W (the closed surge tank transfers heat to the sandstone rock mass) has been observed in the air domain, as shown in Figure 13b. If the static pressure increases up to 200 kPa inside the surge tank, the air temperature increases up to 363 K. The water heat transfer is fully negligible. The obtained results show that the thermodynamic behavior during transient processes depends on the surge tank size, rock mass material, and period of the mass oscillations.
Energies 2019, 12, x FOR PEER REVIEW 13 of 17 Figure 13 shows the total temperature and the heat flux through the walls of the closed surge tank for both water and air phases. The water temperature remains constant throughout the transient process. However, the air temperature increases when the static pressure increases. A heat transfer of -80 W (the closed surge tank transfers heat to the sandstone rock mass) has been observed in the air domain, as shown in Figure 13b. If the static pressure increases up to 200 kPa inside the surge tank, the air temperature increases up to 363 K. The water heat transfer is fully negligible. The obtained results show that the thermodynamic behavior during transient processes depends on the surge tank size, rock mass material, and period of the mass oscillations.  Table 4 and Figure 14 present a comparative analysis of analytical and numerical results for the different guide vanes closing schemes. Maximum static pressure (in relative terms) and water level fluctuations in the surge tank are indicated for both analytical and CFD simulations. The variation obtained between both methods is roughly 7.9% for pressure values.   Table 4 and Figure 14 present a comparative analysis of analytical and numerical results for the different guide vanes closing schemes. Maximum static pressure (in relative terms) and water level fluctuations in the surge tank are indicated for both analytical and CFD simulations. The variation obtained between both methods is roughly 7.9% for pressure values.

Conclusions
In the current energy transition, where VRE are rapidly growing, ESS, such as PSH, have become extremely important to stabilize and balance the global electricity grids. A hydraulic transient analysis and thermodynamic behavior of an UPSH plant in pumping mode was investigated. A closed surge tank with 50 m in height and 5 m in diameter, excavated in a sandstone rock mass at 450 m depth, has been considered. The model also considers a vent shaft with 1 m in diameter. Analytical and 3D CFD numerical simulations have been carried out while considering different linear closing schemes of the guide vanes (1, 10, and 15 s). The variation in static pressure in the air duct, tunnel, lower reservoir, and surge tank, as well as the heat transfer in the closed surge tank, were presented in the transient process.
During the operation in pumping mode, the static pressure in the air duct reaches maximum underpressures (-1.6 kPa). Additionally, the air velocity and static pressure in the air duct vary during the transient process. At the beginning of the closing process, the air pressure in the air duct decreases. The maximum values of the static pressure in surge have been 112.8 kPa for a closure time of 1 s. For the same closing scheme, the maximum water level fluctuation reaches 6.5 m. As observed in the simulations, when the diameter and the height of the surge tank are reduced, the initial volume of air is reduced and, therefore, the pressure and the level of water are also reduced, extending the duration of the transient process. The air temperature in the surge tank increases when the static pressure increases. As verified in the numerical simulations, a heat flux through the surge tank walls of -80 W has been obtained in the existing air in the closed surge tank. When the static pressure in the air phase increases, the temperature and, therefore, the heat flux also increase, varying the thermodynamic behavior. This study demonstrates the paramount importance of the dimensioning of the underground infrastructure of UPSH plants. The obtained results show that the thermodynamic behavior during transient processes depends on the surge tank size, rock mass material, and period of the mass oscillations.
The use of closed mines as underground reservoirs of UPSH plants has important environmental benefits. These initiatives also aid in ensuring economic development of depressed mining areas after closure. However, when compared to conventional systems, there are several drawbacks, such as higher investment and operation and maintenance costs. In addition, the

Conclusions
In the current energy transition, where VRE are rapidly growing, ESS, such as PSH, have become extremely important to stabilize and balance the global electricity grids. A hydraulic transient analysis and thermodynamic behavior of an UPSH plant in pumping mode was investigated. A closed surge tank with 50 m in height and 5 m in diameter, excavated in a sandstone rock mass at 450 m depth, has been considered. The model also considers a vent shaft with 1 m in diameter. Analytical and 3D CFD numerical simulations have been carried out while considering different linear closing schemes of the guide vanes (1, 10, and 15 s). The variation in static pressure in the air duct, tunnel, lower reservoir, and surge tank, as well as the heat transfer in the closed surge tank, were presented in the transient process.
During the operation in pumping mode, the static pressure in the air duct reaches maximum underpressures (−1.6 kPa). Additionally, the air velocity and static pressure in the air duct vary during the transient process. At the beginning of the closing process, the air pressure in the air duct decreases. The maximum values of the static pressure in surge have been 112.8 kPa for a closure time of 1 s. For the same closing scheme, the maximum water level fluctuation reaches 6.5 m. As observed in the simulations, when the diameter and the height of the surge tank are reduced, the initial volume of air is reduced and, therefore, the pressure and the level of water are also reduced, extending the duration of the transient process. The air temperature in the surge tank increases when the static pressure increases. As verified in the numerical simulations, a heat flux through the surge tank walls of −80 W has been obtained in the existing air in the closed surge tank. When the static pressure in the air phase increases, the temperature and, therefore, the heat flux also increase, varying the thermodynamic behavior. This study demonstrates the paramount importance of the dimensioning of the underground infrastructure of UPSH plants. The obtained results show that the thermodynamic behavior during transient processes depends on the surge tank size, rock mass material, and period of the mass oscillations.
The use of closed mines as underground reservoirs of UPSH plants has important environmental benefits. These initiatives also aid in ensuring economic development of depressed mining areas after closure. However, when compared to conventional systems, there are several drawbacks, such as higher investment and operation and maintenance costs. In addition, the possible impacts of droughts and changes in mine water chemical parameters could also affect the operation of this subsurface energy storage systems