TES-PD: A Fast and Reliable Numerical Model to Predict the Performance of Thermal Reservoir for Electricity Energy Storage Units

: The spread of renewable resources, such as wind and solar, is one of the main drivers to move from a fossil-based to a renewable-based power generation system. However, wind and solar production are difﬁcult to predict; hence, to avoid a mismatch between electricity supply and demand, there is a need for energy storage units. To this end, new storage concepts have been proposed, and one of the most promising is to store electricity in the form of heat in a Thermal Energy Storage reservoir. However, in Thermal Energy Storage based systems, the critical component is the storage tank and, in particular, its mathematical model as this plays a crucial role in the storage unit performance estimation. Although the literature presents three modelling approaches, each of them differs in the considered parameters and in the method of modelling the ﬂuid and the solid properties. Therefore, there is a need to clarify the model differences and the parameter inﬂuences on plant performance as well as to develop a more complete model. For this purpose, the present work ﬁrst aim is to compare the models available in the literature to identify their strengths and weaknesses. Then, considering that the models’ comparison showed the importance of adopting temperature-dependent ﬂuid and storage material properties to better predict the system performance, the authors developed a new and more detailed model, named TES-PD, which works with time and space variable ﬂuid and solid properties. In addition, the authors included the tank heat losses and the solid effective thermal conductivity to improve the model accuracy. Based on the comparisons between the TES-PD model and the ones available in the literature, the proposal can better predict the ﬁrst cycle charging time, as it avoids a 4% underestimation. This model also avoids overestimation of the delivery time, delivered energy, mean generated power and plant round-trip efﬁciency. Therefore, the results underline that a differential and time-accurate model, like the TES-PD, even if one-dimensional, allows a fast and effective prediction of the performance of both the tank and the storage plant. This is essential information for the preliminary design of innovative large-scale storage units operating with thermal storage. and F.D.V.; Formal analysis, A.B., F.D.V., E.G., A.S. and G.C.; Methodology, A.B., F.D.V., E.G., A.S. and G.C.; Validation, A.B., F.D.V., E.G., A.S. and G.C.; Writing—original draft, A.B. and F.D.V.; and editing, A.B., F.D.V., E.G., A.S. read


Introduction
The ever-increasing energy demand and the need to replace approximately 40% of the currently in-operation power plants due to ageing means that, by 2040, 7200 GW of new capacity needs to be installed [1]. Since the power industry is responsible for approximately 25% of the global CO 2 emissions [2], and considering that there is the need of knocking down the energy sector greenhouse gas emissions, the vast majority of the new capacity will be, on the one hand, of renewable origin (mainly wind and solar photovoltaic units) and, on the other hand, from low emissions fossil-based thermal units (e.g., combined-cycle units) and nuclear power plants.
However, despite the technology improvements and researchers' efforts, nuclear and combined-cycle power plants lack in flexibility (see, e.g., [3][4][5][6][7][8][9]); therefore, the future power production systems will be made by Variable Renewable Energy Sources (VRESs) units and low flexibility thermal power plants, a generation mix that sets out the installation of Large-Scale Energy Storage Systems. Only spreading these large-scale energy storage installations can help balance the production of VRESs without encountering strong power fluctuations that can result in management and control problems, devices faults and, in the worst scenario, in local or even global blackout.
The International Energy Agency estimates that the United States, Europe, China and India require at least 310 GW of new storage capacity to guarantee the grid stability with a high VRES penetration [10]. In addition, the same international institute also underlined that available storage technologies, such as Pumped Hydro Storage (PHS), Compressed Air Energy Storage (CAES) and batteries, would not be able to cope with the need for storage capacity due to geographical constraints (PHS and CAES) or reduced lifespan (batteries) [11].
Therefore, there is the pressing need to make available storage units able to (i) boost VRES installations, (ii) be installed near both VRES and low flexibility thermoelectric units to reduce transmission losses and network congestion, (iii) ensure a grid safe operation and (iv) do not increase the grid overcapacity. The latter is a fundamental point because, at the time of writing, fossil-based thermal plants are used as backup units to compensate VRES power fluctuations.
Indeed, fossil capacity is not directly replaced by renewable plants, resulting in an installed overcapacity corresponding to a rapid decline of fossil units operating hours [12]. Moreover, additional overcapacity and a further reduction of fossil-fuelled plants working hours encounter in the case of extensive exploitation of large-scale energy storage. Thus, today's better option to overcome these issues is to look beyond the traditional notion of storage systems and switch to the "Integrated Energy Storage System" (I-ESS) concept: storage units embeddable into existing fossil thermal power plants [13,14].
In a nutshell, devices, such as gas turbines, electric generators, step-up transformers and transmission lines, which are commonly mounted on existing thermal power units, can be used to build up the I-ESS plant, ensuring that the storage unit (i) does not add overcapacity, (ii) helps in the revamping of underutilised units, (iii) reduces soil consumption, (iv) boosts VRES spread and (v) guarantees safe network operation even with significant mismatches between electricity supply and demand.
For this purpose, the authors developed an I-ESS plant able to be integrated into underutilised thermal power plants and built near the VRES facilities or in an ad hoc site [13,14]. The plant stores electricity as sensible heat in a high-temperature artificial tank consisting of a packed bed and usually called "Thermal Energy Storage" (TES). The deployment of sensible TES guarantees greater cost-effectiveness compared to latent and thermochemical heat storage systems, as demonstrated by several works published in the literature (see, e.g., [15][16][17][18]).
Compared to PHS, which is the most widespread technology with over 96% of the total installed capacity [19], and CAES, such technology does not suffer from geographical constraints and does not require stable water (like PHS) or natural gas (as CAES) flow rate. A long-cycle life also characterises it compared to batteries with a higher simplicity than Pumped Thermal Energy Storage (PTES) [20][21][22][23][24][25][26][27][28][29]. The I-ESS plant proposed by Benato and Stoppato [13,14] is an open-cycle adopting air as working fluid in both charging and delivering mode.
In a nutshell, the charging scheme is made up of a high-temperature tank, a fan, an electric heater, an electric motor and a heat exchanger, while the delivery unit consists of a gas turbine in which the high-temperature tank replaces the combustion chamber. Compared to the available PTES layouts, the adoption of a unique hot tank instead of two (one at high and one at low temperature) and an open-cycle, which allows the use of air instead of argon, reduces the devices' number, the scheme complexity, and the system management.
Despite a large number of published works and patents concerning TES-based storage units, the vast majority of them aim at the development and modelling of innovative plant configurations for storing electricity as sensible heat, while only a few are focused on tank behaviour investigation. The tank arranged as a packed bed is the core of the TES-based storage system, and its mathematical modelling plays a crucial role in the plant performance prediction as well as in the computation efforts and simulation time. In addition, the non-stationary behaviour of such storage remains a central issue also in its management. In particular, when the TES-based energy storage unit is coupled with VRES plants, it is crucial to accurately predict the charging and discharging times, the stored energy, and the energy that can be delivered, based on the state of charge. These parameters depend not only on the plant layout but also mainly on the physical nature of the storage material and the fluid evolution in the storage itself. Therefore, the availability of accurate storage mathematical models remains a pivotal point for predicting plant performance and its real-world implementation.
For these reasons, the storage mathematical model needs to consider the essential variables required to describe its behaviour accurately without significantly increasing the system complexity and computational time. These needs arise from two facts: (i) the tank model is inserted into the storage unit; therefore, the entire model has to provide accurate results in a reasonable time, (ii) computational Fluid Dynamics (CFD) analysis of the sole tank can provide accurate data, deriving from the solution of the Navier-Stokes system of equations using some closure logic for turbulent phenomena.
Such a strategy would take charge of the local variability of many time-related events and phenomena. However, to be truly effective and provide a complete description of the system phenomena, CFD-based approaches require a three-dimensional tank model: a design condition that clashes with the engineering needs of a reliable, accurate and, to some extend, fast model usable during the design phase.
Note that, due to an accurate and fast simulation of the entire plant, the charging and delivery times as well as the stored and deliverable energy can be estimated in light of storage requirements and packed bed constituting materials. Subsequently, using the CFD techniques, the storage tank can be further investigated to explore the local variability of a large number of time-related events and phenomena.
For these reasons, the work initially focuses on comparing the literature most adopted one-dimensional mathematical models concerning the influence of the single parameters on the tank and storage system performance. According to the authors' best knowledge, three models have been used in the literature. Apparently, all of them can accurately describe the time-varying storage tank behaviour, and they have not previously been compared.
Thus, as a first point of novelty, the different tank models are inserted in the energy storage system architecture presented by Benato and Stoppato [13,14] and compared to predict the performance of such systems. As the models' structures are quite different, our investigation was conducted with the logic of increasing modelling complexity to ensure its effectiveness. At first, the behaviour of the storage tank was analysed through the purely algebraic description proposed by Howell et al. [30]. This modelling strategy, due to its ease of deployment, is one of the top choices in the literature, but it is still the most limited.
The second modelling approach is based on the proposal made by McTigue et al. [31]. The model feels a pseudo-time-varying approach, for which the non-stationary of the storage system is treated in terms of average quantities, i.e., the mean temperature and average mass-flow rates. A partial variability of both the operating fluid and the storage solid's thermal coefficients, albeit with a time interval averaged logic, is taken into account. The last modelling approach is the one proposed by Desrues et al. [21]. The model describes the reservoirs' transient behaviour through a purely one-dimensional time marching partial differential equations approach.
The comparison among the three models allows us to underline the importance of (i) space and time variability of the thermophysical coefficients related to the operating fluid and the storage material, (ii) the overall heat loss coefficient and (iii) the solid effective thermal conductivity. Thus, starting from the formulation proposed by Desrues et al. [21], the authors developed a new tank model (named TES-PD model), which includes all these features.
Despite its one-dimensionality but due to its high level of completeness and ability to take charge of space and time variability of the fluid and storage material thermophysical properties, the TES-PD tank model is of real help in designing and managing such devices by predicting data that can be otherwise unobtainable; the latter can then be used as preliminary designing guidelines or input for a subsequent CFD simulations. In addition, the completeness of the proposed model makes it possible to investigate storage systems under innovative operating conditions, such as the one at low temperatures, a research field still partially unexplored.
The work is organized as follows: in Section 2 the TES-PD model is discussed in terms of its analytical description and numerical treatment. Section 3 presents the validation procedure of the models available in the literature in order to use them to highlight the peculiarities of the TES-PD model. Section 4 describes the plant scheme used to test the proposal, while the investigation outcomes and the comparison between the literatureavailable models, and the TES-PD one are summarised and discussed in Section 5. Finally, Section 6 states the conclusions.

The TES-PD Numerical Model
In this section, the TES-PD tank mathematical model used to predict the storage behaviour is presented. The model rearranges the most used tank models and improves their predictive capabilities with additional features. Those models, in order of complexity, are the Mumma & Marvin model [30], the Schumann & McTigue et al. model [31] and the Desrueset al. model [21]. For the sake of conciseness, such numerical strategies are reported in Appendix A and, in the following, they are addressed as "model m.1", "model m.2" and "model m.3", respectively.
The TES-PD model relies on the formulation proposed by Desrues et al. [21]. Despite the common starting point, the authors improved it by the integration in the TES-PD model of the entire set of parameters that affect the behaviour of the storage tank. These parameters are variously missing in the available literature contributions; in particular,

1.
The space and time variability of all the thermophysical properties of both fluid and storage material.

2.
The overall heat loss coefficient.

3.
The solid material effective thermal conductivity.
Model m.3 was selected as the basis for the TES-PD model since it represents the most advanced description of the storage tank already available in the literature. However, model m.3 envisages constant solid and fluid properties computed at an arbitrary reference temperature. This limitation often leads to mispredicting the system dynamics and introducing uncertainties related to the reference state. Moreover, it can not fit the usage of the most advanced storage materials.
If adequately modified, model m.3 allows updating all the thermophysical properties in both space and time. Thus, in the TES-PD model, the authors embedded the update of both the solid and the fluid properties with a procedure similar to the one described by McTigue et al. [31] for model m.2. However, in contrast to model m.2 for which no specifications were reported as far the computation of the thermophysical coefficients, solid material properties were acquired from the National Institute of Standards and Technology database [32], while the fluid thermophysical properties were taken from the CoolProp library [33]. In this manner, the thermophysical properties variability was considered layer-by-layer, updating their values during the time-marching procedure as a function of the actual temperature.
Moreover, concerning the overall heat loss coefficient, this is accounted for only by model m.1 (see Appendix A.1 for the details). However, such a parameter is of fundamental importance for a storage tank because the charging and discharging periods can last from hours/days, making this a non-trivial contribution concerning its performance. Introducing this hypothesis, the TES-PD model reads as follows: The above-mentioned set of equations consists of the fluid's mass conservation equation, the fluid energy conservation equation and the solid energy conservation equation. The model is closed by a constitutive expression for the pressure drop and the ideal gas equation of state: where t is the time and x is the axial tank coordinate with L its length; ρ f , v f , T f and p are the fluid density, velocity, temperature and thermodynamic pressure, respectively; while ρ s and T s are the storage material density and temperature, respectively. k f and k s,e f f denote the fluid and the solid effective thermal conductivity, and c p, f , and c p,s are the fluid and the solid specific heat coefficient at constant pressure, respectively. ε denotes the void fraction, T amb is the ambient temperature, and r is the specific gas constant. The α and β coefficients generalise model m.3, assuming different values depending on the storage material geometry. α equals to 6d −2 Nu k f c −1 p, f for spheres packing and equals to 4εd −2 h Nu k f c −1 p, f for microchannels. d denotes the particle equivalent diameter, while d h denotes their hydraulic diameter. As far as the β coefficient is concerned, the latter equals (1 − ε)(dε) −1 for spheres packing, while it equals 4d −1 h for microchannels. C u (1 − ε) −1 denotes the external surface to volume ratio, while U i is the overall heat loss coefficient.
C u values are 4L −1 and 2π A −1 √ Aπ −1 when the considered geometry is a microchannel structure or sphere packing, respectively. Finally, the hydraulic diameter, d h , together with the Nusselt number, Nu, and the friction coefficient, C f , are computed according to Desrues et al. [21] as given in Appendix A.3. In the present work, the authors considered packed sphere geometries where stated.
Compared to model m.3, Equation (1a,b) remain formally unaltered, while Equation (1c) accounts for the time-variability of the thermophysical coefficients and embeds the heat losses through the external surface. Such improvements are of fundamental importance in detailing the system's dynamics since truly layer-by-layer variability of the solid and fluid properties can be accounted for as well as the heat exchange between the tank and the external ambient temperature.
To solve the previously-mentioned set of partial differential equations, the authors adopted a first-order explicit Euler scheme for the temporal integration, while a firstorder generalised upwind finite difference method was employed to compute the first spatial derivatives. A second-order centred finite difference approach was used for second derivatives treatment [34,35]. The cell-bound variability of the physical coefficients was accounted for through the method proposed by De Vanna et al. [36].
The fluid velocity was updated through a segregated approach to overcome the equation coupling. This set of discretization schemes is guided by the simplicity and computational efficiency. While using first-order methods, this choice proved to produce accurate results, which fits the available references in corresponding conditions. Note that the literature does not suggest methods to solve the set of equations; hence, the authors developed the above-mentioned approach and integrated it into the TES-PD model in order to provide both an updated TES model and a simple and computationally efficient resolution scheme.
Concerning the initial conditions for the partial differential equations problem, the entire set of tests was performed by the authors while adopting uniform ambient temperature and density, and the flow was considered initially at rest. The fluid properties were set according to the ideal gas law, while the values for the fluid mass flow rate, pressure, and temperature were enforced at the inlet location depending on the plant cycle in which the tank was installed.
Zero-gradient conditions, with first-order of accuracy, were set for the solid's temperature distribution at the tank's inlet and outlet. The strategy allows modelling the storage material as adiabatic, packing the heat losses in the dispersion coefficient only.
The modelling hypothesis introduced in the TES-PD model was fully validated. In particular, the proposed strategy produced similar heat losses to the one predicted by model m.1 (if used in corresponding conditions), the latter being the only model that accounts for this contribution. Several analyses were also performed looking at the role of non-constant thermophysical properties, tracing their contribution to the overall tank performance. More details of such investigations are given in Sections 3 and 5.

Models Validation
Before using the TES-PD model to predict the storage tank dynamics and, subsequently, the I-ESS performance, it is fundamental to validate it. However, to the authors' best knowledge, the literature lacks experimental data as well as numerical results obtained with detailed models like the proposal. For these reasons, the models m.1, m.2 and m.3, which partially enclose some of the TES-PD model features, can be used to demonstrate the accuracy and computational efficiency introduced by proposed model.
To this end, there is a need to implement those models and validate them in order to ensure the quality of the further comparisons. In particular, Section 3.1 describes the validation process of the model m.2 while Section 3.2 presents the procedure adopted to validate the model m.3. The validation process of model m.1 is not discussed in this work since it has been extensively analysed in [13,26], while the mathematical description of model m.1, m.2 and m.3 is summarised in Appendix A.

Validation of Model m.2
To validate the reservoir model according to the hypothesis of McTigue et al. [31], there is a need to consider their PTES plant arrangement [37]. The latter comprises one pair of turbomachines, two heat exchangers, and two packed beds: one at high pressure and temperature (named "hot tank") and one at ambient pressure but low temperature (called "cold tank").
During the charging phase, the system acts as a high-temperature heat pump cycle where the electricity is first converted into heat through the compressor and then stored as heat and cold in the hot and cold reservoirs, respectively. During the discharge, the plant acts as a closed Brayton cycle heat engine; thus, the system allows for extracting the stored energy and re-generating electricity. Both the hot and cold tanks are characterised by a height L of 5 m and a diameter D of 5 m; thus, L/D = 1.
The void fraction, ε, is set equal to 0.33, while the particle diameter d is assumed to be 0.03 m. The working fluid is argon and the mass-flow-rate,ṁ, is set equal to 12.5 kg·s −1 . The pressure ratio between the two tanks is set equal to 10. Further details of the plant description and modelling assumptions can be found in [31]. Figure 1 shows the temperature trends inside the cold tank, considering a stabilized cyclic operation. The charging and discharging periods depend on a crucial parameter for cyclic operations: the utilization factor Π. This parameter, equal to 0.75 for both the phases, indicates the ratio between the charging and nameplate time (based on the thermal front speed estimation in ideal conditions).
Temperature trends refer to three different operating conditions during the charging process, considering constant or variable the solid storage material's specific heat c ps . As given in Figure 1, the depicted trends are in good agreement with the reference for both constant and variable solid properties. In particular, the initial stages' temperature trends show a maximum difference of 5%, while, at the end of the process, this decreases to almost 1%. The uncertainty for some settings (e.g., the update frequency or the component models) and material property values may affect the results, and the slight mismatches are to be ascribed to these factors.   . Temperature distribution inside the cold reservoir during the charging process, at different charging stages (5%, 50% and 95%). Our are compared with the results reported in [31,37]. In black empty dots, temperature profiles are reported considering constant the storage material's specific heat (the mean value is selected). In black filled dots, temperature profiles concern the case of a variable storage material's specific heat: each layer has a different value for this parameter, according to the actual temperature.
To further validate the implemented model, another possible configuration for this PTES system emerges from the optimisation study performed by McTigue et al. [31]. In the reference, a more realistic scenario arises since the design implies the use of a polytropic efficiency equal to 0.9. The graphical extraction of the geometrical properties leads to a L/D ratio of 0.57 and 0.23 for the hot and cold reservoir, respectively, a void fraction ε of 0.35, and a particle diameter of 0.63 cm for both the tanks. The pressure ratio during the charging phase was maintained equal to 10.92, while during the delivery, it decreased to 9.105. Figure 2 shows the round-trip efficiency and the energy density (i.e., the ratio between the restored energy and the total volume of the reservoirs) trends for the first 40 cycles, employing a utilisation factor of 0.25. The steady-state is achievable after many more cycles due to the lower utilisation factor Π. The reported trends for these two parameters consider a constant as well as a variable solid's specific heat. The comparison involves the values reported by McTigue et al. [31], which were 0.59 for the efficiency and 65 MJ m −3 for the energy density.
Since McTigue et al. [31] sized the PTES system to produce an ideal energy density of 300 MJ m −3 ; the results in terms of energy density are in full agreement with the available data.
Hence, the model m.2 can be considered fully reproduced because the temperature profiles, the round-trip efficiency and the energy density trends were replicated concerning boththe constant and variable solid specific heat at constant pressure, c p,s .  [31]. (a) The round-trip efficiency obtained by the authors is compared with the one reported in the reference. (b) Energy density, defined as the ratio between the energy re-obtained during the discharging phase and the tanks' total volume. The latter agrees with that reported in the reference since the energy density is equal to 300 MJ m −3 in design conditions.

Model m.3 Validation Process
To validate the bed model according to the model m.3 description, it was necessary to build up the PTES system proposed by Desrues et al. [21]. The latter is similar to the McTigue et al. [31] one, with some changes in the components' management and cycle design parameters.
The numerical setup, in terms of the grid and time-step, accounts for the following choices: two hundred layers were used for the spatial discretisation allowing for a spatial resolution of ∆x = 0.05 m. Such resolution was found to be sufficient in guaranteeing grid independence, and the reader is directed to Appendix B for further details. Temporal stability is recovered setting the time step, ∆t, equal to 0.01 s. The storage material is considered made of microchannels; thus, the α and β coefficients were selected, consequently. The global cycle update time (i.e., all the inputs and outputs update of each component of the system) was set equal to 120 s.
The process' steady state was reached after approximately ten cycles, considering a cyclic operation of charge-discharge. The pseudo-stationary results were compared to those reported in the reference (graphically extracted). The plant was managed by controlling the temperature at the outlet of the tanks; thus, the utilisation factor was not defined. Figure 3 shows the trends of stored and delivered energy. Stored energy represents the whole energy absorbed by the system during the charging phase, while delivered energy is the net energy output during discharge. Another trend shown is the efficiency over cycles (i.e., the ratio between restored and stored energy), which is in good agreement with the reference values.
The mismatch in energy between the authors' implementation and the reference [21] depends on the not explicitly reported choices: for example, the storage material used or the other component model. Thus, the efficiency trend is more representative of the numerical model adopted since it is more independent of such particular choices. Cycles Stored and restored energy  Figure 3. Energy and efficiency during the first 10 charge-discharge cycles for Desrues [21] and the present case. (a) In black are reported the authors' stored energy (empty dots) during the charging process and the delivered energy (filled dots) during the discharging process. In red are shown the results of the reference. (b) Dots in black indicate the efficiency in the authors' case, defined as the ratio between restored and stored energy. These trends show good accuracy in reproducing the process efficiency, compared with those described in the reference (red empty squares). Figure 4 also illustrates thermal fronts inside the hot and cold tanks. The reported results refer to the 10th cycle at the end of both the charge and discharge phases. As clearly depicted, even these trends show a good match concerning the reference solution, particularly where the temperature gradient is high. Thus, the authors convey that the model m.3 implementation can be considered validated since the maximum temperature difference occurs at the tank limits (mostly influenced by other components' simulation choices), while the internal area gradient shows a less than 1% difference.
Considering that the TES-PD model is based on model m.3, the authors introduced their model in the plant scheme proposed by Desrues et al. [21], and, by avoiding the use of variable thermophysical properties and heat losses through the surface of the tank, they verified that the obtained results were exactly the ones obtained with model m.3. Therefore, it is realistic to consider the TES-PD model validated under model m.3 conditions. In any case, to support this claim, further investigations will be presented in Section 5 with regard to the selected test case.

The Selected Test Case
To carry out an objective comparison among numerical reservoir models, there is a need to include each of them in a plant. For this purpose, the I-ESS configuration developed by Benato and Stoppato [13] has been selected. The plant model is linked to the National Institute of Standards and Technology [32] database and CoolProp [33] library to acquire the storage material and the working fluid thermophysical properties. More information on plant component modelling and validation procedure were reported by Benato and Stoppato [13] and Benato and Stoppato [14], while, in the following, a brief description of the plant scheme and control strategy are given.
The plant layout used to store electricity in the form of thermal energy is depicted in Figure 5a, while Figure 5b illustrates the discharge configuration that allows converting the thermal energy into electricity.
The charging layout consists of a fan (FAN), a heat exchanger (HX), an electric heater (EH) and a high-temperature storage tank as shown in Figure 5a. The plant acts as an open-cycle where the air is used as the working fluid. The fan sucks the air from the ambient atmosphere and forces it to follow the path from 1 to 6. The electricity generated in surplus, e.g., by the VRES plants, is converted into heat through the electric heater.
In practice, the device uses electrical energy to increase the air temperature from T 3 to T 4 . The latter is equal to the maximum cycle temperature, T max , and it is a design parameter. Compared to other TES-based layouts where the electricity conversion into heat is done by the compressor, the electric heater allows maintaining T 4 constant and equal to the maximum one independently to the value assumed by T 3 .
After the conversion, the electricity is stored as thermal energy in the hot storage, which made by a packed bed, is initially at ambient temperature and, during the charge, is slowly heated up by the hot air entering it at a temperature equal to T max . Since the plant is based on an open-cycle, a heat exchanger is inserted to reduce the heat losses. The device, if needed, should be bypassed and installed only in storage units characterised by medium-to-low size to contain its dimensions and, subsequently, the costs. Note that the charge stops when the mismatch between the maximum temperature, T max , and the air temperature at the tank's outlet, T 5 , is equal to a pre-computed value, named the charging tolerance (tol). The charging tolerance (which can be equal to, e.g., 10, 100 or 300 K) widely affects the performance of the plant in terms of the charging and delivering times (t CHG and t DIS ), stored and delivered energy (E CHG and E DIS ) and round-trip efficiency (η rt ). Therefore, based on the service requested to the I-ESS plant, the charging tolerance needs to be carefully selected and, in the case of switching the operation mode, reassessed. Figure 5b depicts the delivery arrangement and shows the route follows by the air: from 11 to 16. The system is composed of a compressor, an expander, an electric motor/generator (MG), a heat exchanger (HX) and a high-temperature tank. The generation unit acts as a modified gas turbine in which the combustion chamber has been replaced with the high-temperature storage tank.
As previously stated, the idea of configuring the I-ESS plant using gas turbine components stems from the fact that gas turbines operating in open-cycle are not able to compete in the actual electricity market due to high operating costs and low efficiency. However, in this scenario, their utilisation hours dropped rather steeply. Thus, gas turbines conversion in storage units can be helpful to revamp such plants as well as to redeploy industrial sites.
The compressor sucks the air at ambient conditions and increases its pressure and temperature. Then, the fluid passes through the heat exchanger and enters the hot tank, where it is heated up. At the beginning of the delivery, the air leaves the hot tank at a temperature approximately equal to T 14 = T max − tol. At this point, the air is expanded in the turbine and injected into the heat exchanger. The expansion process generates mechanical power. Part of it is used to drive the compressor while the remaining one is converted into electricity through the electric generator.
After the expansion, the air is cooled down into the heat exchanger and, then it is released into the environment. As for the charging scheme, the heat exchanger should be bypassed if needed and can be installed only in storage units characterised by mediumto-small size to contain its dimensions and, subsequently, costs. As is evident, the heat exchanger and the tank shown in Figure 5a are the same components of Figure 5b.
As for the charge, the delivery phase can also be stopped at different points depending on the service requested to the I-ESS plant. The state of discharge influences the time requested for the second to i-th charge, the energy that can be stored during the second to i-th charge, the discharge time and the round-trip efficiency. In addition, during the discharge, the generated power (P GEN ) by the I-ESS plant decreases with the reduction of the tank temperature. Nevertheless, to guarantee an almost constant generation during the delivery phase, the plant control system stops the discharge when the mismatch between the generated power and the design one is higher than 5%.
The tank's size has to correspond to the amount of energy that needs to be stored at the required temperature. Therefore, as suggested by Singh et al. [38], the size of the packed bed needs to be fixed so that, during the charge, the bed absorbs the maximum amount of energy made available by the heat transfer fluid flowing through the bed.
In this way, at the beginning of the charge, the fluid leaves the storage at a temperature approximately equal to the initial one. Conversely, at the end of the charge, the bed can achieve a temperature near the one of the fluid entering the storage. Note that, during the charge, the heat transfer fluid flows from the top to the bottom, while, during delivery, the flow is reversed. The tank is cylindrically shaped with an upper plenum, a packed bed storage and a lower plenum as depicted in Figure 5, while it is arranged vertically to avoid buoyancy-driven instabilities of the thermal fronts.
We selected an existing tank placed in an underutilized facility that can accommodate a packed bed with a volume of 500 m 3 and a height of 8.30 m. We assumed that the tank temperature at the beginning of the first charge is equal to ambient temperature, 25°C. The packed bed void fraction was assumed equal to 0.4 because the tank is filled with spheres of aluminium oxide (0.05 m in diameter) with a density and a mean specific heat equal to 3990 kg m −3 and 1150 J kg −1 K −1 , respectively [39].
The authors selected the above-mentioned packed bed specifications following the indications given by Ataer [40] and Singh et al. [41]. As suggested by Desrues [42], for this geometrical arrangement, the solid material's effective thermal conductivity can be assumed up to 1 W m −1 K −1 , while the overall heat loss coefficient U i can be set equal to 0.7 W m −2 K −1 .
Considering again the idea of converting underutilized thermal units into a storage system, the authors, with the help of an Italian gas turbine manufacturer, selected a 3 MW gas turbine unit installed in an Italian industrial site as a device to be used for the delivery phase. Its operating parameters included a mass-flow rate of 16 kg s −1 , a pressure ratio of 9.5 and a maximum temperature at the turbine inlet of 1123 K. Therefore, the maximum temperature of the cycle, T max , was fixed equal to 1123 K.

Results and Discussion
Before proceeding with the model comparisons, it is crucial to remark that each model is not different from the others concerning the structure only; however, there are even different modelling choices. Model m.1 considers the heat losses, while model m.3 considers the solid thermal conductivity contribution but it is not reciprocal. Model m.2 does not consider any of these effects. In addition to these aspects, model m.1 and m.3 required constant properties for both solid and fluid, i.e., c p,s , c p, f , k f and µ f while model m.2 and model TES-PD introduce the property variability.
Intending to understand the influence on the plant performance of both the heat loss coefficient (considered in m.1) and the effective thermal conductivity (considered in m.3 and TES-PD), the authors first compare the models overriding these contributions. However, considering that models m.2 and TES-PD required variable properties and that model TES-PD becomes equal to m.3 when the constant properties are accounted for and the heat loss coefficient and the effective conductivity neglected, the authors ran models m.1, m.2 and m.3 examining different temperatures at which the fluid and the solid material properties were computed.
Note that the authors imposed the same constant property behaviour even for model m.2, aiming to compare the results. After that, the obtained findings can also be compared with the outcomes retrieved from model m.2 and TES-PD model run with variable properties of the fluid and the solid. The outcomes of this analysis are given in Section 5.1.
In Section 5.2, instead, we investigated the role of the charging tolerance in the tank performance. The dynamic performance of a TES-based storage plant was affected not only by the fluid and solid properties but also by the heat losses and effective conductivity and charging tolerance. The lower the mismatch between the maximum temperature and the air temperature at the tank's outlet was, the higher the charging time and the stored energy were.
Therefore, it is interesting to compare the models' behaviour also with different charging tolerance. To overcome issues related to fluid and solid properties evaluations, the comparison among the models was made considering constant properties also in the case of model m.2 and neglecting the heat loss coefficient and the effective conductivity. Under these conditions, model TES-PD falls back to the results obtained from model m.3.
Finally, in Section 5.3, the models are compared considering their original formulation to underline the pros and cons of a more complete formulation. Regarding the storage performance affected by the solid material conductivity and given the lack of experimental data, the authors performed a sensitivity analysis based on literature data and available correlations.
As far as the simulations setup, for each model, a specific layer height was based on the solution sensibility considerations, ensuring less than a 2% difference in the observed results, doubling the number of layers. Thus, the chosen values of grid refinement were ∆x equal to 0.1 m (83 layers, [13]) and 0.0166 m (500 layers, [37]) for model m.1 and m. Comparisons were made analysing the plant performance parameters: charging and delivering times, stored and delivered energy, generated power and round-trip efficiency evaluated for the first cycle, with the first stable cycle and the average cycle calculated by averaging 50 stable cycles.

Influence on Fluid and Solid Properties of the Reference Temperature
Dealing with constant properties leads to assuming a temperature corresponding to the properties' reference value. Admissible reference temperatures can be the ambient temperature, the maximum temperature of the cycle or a mean value between them. In this analysis, due to the absence of literature specifications on reference temperature selection, the authors tested three values to underline the differences introduced by this assumption. Table 1 lists the results obtained when the fluid and storage material properties c p,s , c p, f , k f and µ f were computed with a reference temperature equal to T amb = 25 • C, T max = 850 • C and T mean = 437.5 • C. We also list the outcomes of simulations computed using properties variable in space and time for the models that allow these features (m.2 and TES-PD).
T amb columns list the system performance parameters when setting the reference temperature for properties evaluation to the ambient temperature. The maximum difference involves the first cycle's round-trip efficiency, whose relative difference is about 7% considering different models. Then, in a stabilised cycle, the discrepancies become more relevant for discharging times and energy. T max columns report the system performance parameters when the reference temperature for properties evaluation was set to 850 • C.
The maximum differences drop from 7% to 4% concerning the first charge, compared to the T amb case. In a stabilised cycle, the differences decrease to almost 3.5%, which is a minor mismatch. When assuming T = T mean = 437.5 • C as the reference temperature, the relative differences are almost comparable to the T = T max reference temperature case. Table 1. Comparison among models when the fluid and solid properties (c p,s , c p, f , k f and µ f ) are computed at a fixed reference temperature (T amb , T max and T mean ) and at a temperature, T var , reliant on the layer and time (only for m.2 and TES-PD). In order to compare the models, the heat loss coefficient U i and the solid effective thermal conductivity, k s,e f f , were both set equal to 0, while the charging tolerance was set equal to 10 K. The plant performance parameters are listed for the first cycle, the first stable cycle and average one. Observing these simulation outcomes, it is evident that the differences among the models at a specific reference temperature are negligible. Model m.3 always predicted a higher efficiency, even if the differences were generally variable but, still, less than 9%. This concerns only the differences among the models. More meaningful is the evaluation of how the choice of the reference temperature affects the results.
In this case, since models m.2 and TES-PD can deal with variable properties, it is also interesting to add these results to the discussion (columns T var ). Note that, assuming variable properties means not only updating values during the time-marching procedure but also layer by layer. The solid and fluid's specific heat and the air's viscosity and conductivity were evaluated on the single layer's local temperature basis.
From the obtained results, it is clear that the selection of the reference temperature is a fundamental step because adopting T amb instead of T mean or T max caused an underestimation of the first and average charge stored energy (E CHG ) of 37% and 42%, respectively. Similarly, the energy deliverable during the discharges (E DIS ) was underestimated by a value higher than 47.5%.
Underling these underestimations in the energy storage capacity are essential since it predominantly affects the design volume of the storage reservoir. Independently to the adopted model, selecting T amb as reference temperature instead of T mean or T max can lead to the selection of a storage volume 50% smaller than the one requested from, e.g., the VRES plant. In addition to the above-mentioned design issue, the reference temperature selection also affects the time required to charge, t CHG and discharge the storage, t DIS .
In fact, the adoption of T amb as a reference reduced the first and the subsequent charges times of 33% and 37%, respectively, while the discharge time was underestimated by 47.7%. These underestimations of the charging and the delivery times are the source of wrong evaluations of storage plant market availability and response time: parameters that need to be accurately estimated to avoid storage plant inability of compensating the mismatch between, e.g., the VRES supply and network demand, which can be the cause of devices fault or, in the worst scenario, grid blackout.
Analysing the columns of Table 1, which refer to T var , and comparing them with the outcomes of the simulations derived by the selection of T mean or T max as the reference temperature, it is clear that the mismatches in the charge and discharge times and energy were lower than 3%. Therefore, it appears that adopting a constant temperature equal to T mean or T max for properties computations did not introduce relevant discrepancies.
However, since the choice of any reference temperature is purely arbitrary and, considering that the simulation time with variable properties increases by 10% compared to constant properties at T mean or T max , the use of variable properties can lead to a genuinely consistent model especially in view of including in it parameters, such as heat losses and solid thermal conductivity, as done in the TES-PD model. Figure 6 shows the temperature trends for these tests, adopting model m.3. Since differences between models in this configuration are not too marked, the figure focuses on a single model results' comparison with different reference temperature choices. The most significant difference concerns the results derived using the environmental temperature as the reference.
The charge and discharge processes last almost 30% less than in the other cases, thus, involving less energy. The temperature profiles inside the tank at the end of the first stable cycle's charge and discharge were not greatly affected.
A minor difference concerns the last layer temperature during the first charge, which increased more gradually when the properties were variable. These considerations lead to say that, for these specific conditions, i.e., materials and system behaviour, T max or T mean are good choices as the reference temperature. However, this can not be generalized. Adopting variable properties, a possibility offered by model TES-PD, is a solution that allows avoiding wrong temperature choices, although it is 10% more computationally expensive.

Influence of the Charging Tolerance
The temperature used as the reference for fluid and solid properties computations is not the unique parameter that directly affects the performance of TES-based energy storage systems. In fact, the time requested for charging the storage tank and the energy stored in it is directly linked to the charging tolerance, tol, namely the mismatch between the maximum temperature and the air temperature at the tank's outlet. The lower the tolerance is, the higher the charging time is.
As for the properties' reference temperature, the charging tolerance is also a design parameter whose selection depends exclusively on how the plant is operated, which, in turn, relies on the network or VRES plant needs. For this purpose, we investigated the difference in the system dynamics caused by different the charging tolerance and adopting a tol equal to 10, 100 and 200 K. To better detect model mismatches concerning only the charging tolerance, simulations were performed using model m.1, m.2 and m.3, computing the fluid and solid properties at T mean and neglecting the heat loss coefficient and the solid effective thermal conductivity. Table 2 summarises the simulation outcomes.  Table 2. Comparison between models when assuming different charging tolerances. The charging tolerance, which is the temperature difference between the cycle's maximum temperature and the air temperature at the tank's outlet, was set equal to 10, 100 and 200 K. Even if the maximum relative difference between the characteristic parameters slightly increased (from a 4% to 7% difference in the discharging times and energy during a stabilized operation) when the charging tolerance was 200 K, the most noticeable differences concerned the system's round-trip efficiency, η rt . In particular, the higher the tolerance is, the higher the round-trip efficiency is. In contrast, the lower the tolerance is, the higher the charging and delivery times and the charging and delivery energy are.
Apart from these aspects, comparing the simulations results for a given charging tolerance revealed minor mismatches among them: usually, less than 6% and 8% in the first and first stable cycle parameters, respectively. Therefore, the model comparison does not underline a significant mismatch in the plant performance prediction if the same charging tolerance is set.
However, looking deeper into the results listed in Table 2 reveals that the charging tolerance influenced the stabilisation of the process. When tol equalled 10 K, the performance parameters at the first stabilised cycle were identical to the ones obtained from the average cycle calculated by averaging 50 stabilised cycles; this was no longer valid when setting the charging tolerance to 200 K. The different temperature profile at the end of the charge caused a delay in reaching the stabilisation. The delay effect was more significantly remarkable when the tolerance was higher. In fact, when tol = 10, the temperature was almost T max in the whole tank at the charge's end, making the second discharge process similar to the first one.
To underline this critical aspect, Figure 7 shows the temperature trends over time, adopting the m.3 model with different charging tolerances. Even if the discharge process ends when the generated power is 5% lower than the maximum, the charging tolerance reduction affects this phase. The temperature distributions inside the tank (Figure 7b) show the lower energy availability when setting the tolerance to 200 K. The data in Table 2 state that when tol = 10 K, the average cycle's charging energy was about 190 MWh, compared to around 130 MWh when tol = 200 K.
The same happens to the discharge energy but less markedly since the efficiency slightly increases. In conclusion, the choice of a particular tolerance in the charging phase did not have an exclusive impact on the charging time; however, the whole system's dynamics were conditioned. This aspect is fundamental because such a tank is hooked up to a highly variable energy source, as the system's history heavily affects its overall performance.

Models' Original Formulation
As previously stated, in this work, we aim to compare the available packed bed models, understand the influence of the constituting parameters, and then develop a new mathematical structure to take into account the most impacting parameters while preserving the reliability and computation efficiency. However, to evaluate the impact on plant performance of the temperature at which the fluid and the solid properties are computed (Section 5.1) and of the charging tolerance (Section 5.2), the models need to work under the same operating conditions. Thus, the authors set to zero the heat loss term (U i ) considered in model m.1, the effective thermal conductivity (k s,e f f ) deemed in model m.3 and the property variability through layers in model m.2. In this manner, the analysis revealed that, for models working with variables computed at a constant temperature, a good compromise was to consider T mean (the average temperature between the ambient and the maximum ones) despite the best choice is a model considering layer-by-layer fluid and solid properties' variability (like model m.2 and TES-PD).
This analysis is crucial in light of a model comparison based on their original formulations. In fact, to carry out this analysis, the properties for model m.1 and m.3 were computed considering T mean while m.2 and TES-PD assumed variable properties for each layer, according to the local temperature. The heat loss coefficient U i was set equal to 0.7 Wm −2 K −1 , while k s,e f f was assumed equal to 1 Wm −1 K −1 [30,42]. Regarding the charging tolerance, the authors selected 10 K because they assumed to fully charge the hot storage. Table 3 lists the main outcomes of the performed investigation. Table 3. Test of all models in their original formulation. In this case, the properties for m.1 and m.3 were assumed constant (evaluated at T mean ), while m.2 and TES-PD assumed variable properties for each layer, according to the local temperature. The heat loss coefficient U i was set equal to 0.7 Wm −2 K −1 , while k s,e f f was 1 Wm −1 K −1 . The charging tolerance was fixed equal to 10 K. Compared to values reported in Table 1 (column related to m.1 adopting T mean ), the consideration of the heat loss coefficient increases the first and average cycle charging times of 6.2 and 7%, while the first and average cycle discharging times were reduced to 8%: a fact that guarantees increasing the energy stored in the tank by 2.2% (first cycle) and 1% (average cycle), but that reduces both the discharge energy and round-trip efficiency by approximately 10%. The plant's mean generated power resulted as reduced by approximately 1.2%.
Thus, the heat loss coefficient is a parameter that needs to be taken into account in the tank model to correctly estimate its dynamic behaviour and, subsequently, the storage plant performance. In a nutshell, neglecting U i implies a 10% overestimation of the plant deliverable energy, a fact that can be the source of storage plant management issues.
Compared to results listed in Table 1 (column related to m.3 adopting T mean ), a model formulation in which the solid conductivity was set equal to 1 Wm −1 K −1 as suggested by Desrues [42], did not introduce considerable differences. For the first and the average cycle, the charge and discharge times as well as the charged and delivered energy, the gaps were less than 1%. Therefore, a model formulation in which the solid conductivity was considered constant did not drastically affect the prediction of plant performance.
More interesting is the comparison among the computations performed in Table 1 (column related to TES-PD adopting T var ) and the ones given in Table 3. They reveal that neglecting heat losses and solid effective thermal conductivity largely influenced the plant performance prediction. The first cycle charging time resulted as 4.3% lower while the delivering one was reduced by 6.1%.
The stored and delivered energy were 1.6 and 7.7% lower while the plant round-trip efficiency resulted as 9.4% lower in the first charge-discharge cycle and 8.3% considering an average cycle. Thus, the original formulation of the TES-PD model guarantees considering the account property variability, heat losses, and solid conductivity without compromising the simulation time, which was only 12% higher than for the other models.
For the sake of completeness, in Figure 8, we present the non-dimensional temperature trends in the hot tank during the first and the average charge, the average discharge and the temperature distributions inside the tank at the end of the average cycle's charge (CHG) and discharge (DIS) for the four models.
The trends did not substantially differ among each other, although the completeness of the TES-PD model guaranteed a better estimation of the overall performance. However, in previous simulations, the thermal conductivity value was assumed constant and equal to a value suggested by Desrues [42] despite that solid thermal conductivity can lead to very different results if it is properly modelled. In particular, in contrast to other models, the TES-PD allows the implementation of correlations in which the solid thermal conductivity is variable according to a temperature-based law (see, e.g., [43]).
However, in this case, due to the physical process behaviour, the thermal conductivity can assume a value up to 7, a fact that imposes a drastic reduction of both the simulation time-steps and layer thickness, parameters, which, in turn, affects the computation time, which can be to three orders of magnitude higher. For these reasons and considering the literature lack of accurate data for solid thermal conductivity in the simulated conditions, the authors performed a parametric analysis (see , Table 4), which demonstrated low discrepancies among k s,e f f equal to 0, 1 or 2 W m −1 K −1 .
In contrast, the computational time increased by 20% from 1 to 2 W m −1 K −1 . Therefore, despite the model's ability to take into account the effective conductivity variability according to a temperature-based law, the lack of objective data and the vast increase in computational time led the authors to choose, for the design phase, a constant value equal to 1 W m −1 K −1 . In future TES-PD updated versions and CFD investigations, the effect of this parameter needs to be further investigated to sharpen the model's ability to predict storage dynamic behaviour.

Conclusions
The environmental impact of power generation units is in boosting wind and solar photovoltaic plant production because they can guarantee green electricity generation. However, such plants are characterised by a variable and intermittent nature, a feature that results in considerable power fluctuations, which, in turn, can provoke grid control issues and device faults. To overcome these problems and guarantee a rapid transition from a fossil-based generation mix to a renewable one, installing large-scale energy storage units is required.
For this purpose, we developed a system able to compete with pumped hydro energy storage, compressed air energy storage and batteries. However, as the proposed storage system is based on the thermal energy storage concept, there is a need to predict the storage performance correctly. This is not an easy task considering that the literature reported three tank models, each one different from the others in terms of how the fluid and storage material properties are computed and the consideration or not of the overall heat loss coefficient and the effective solid conductivity. Therefore, to highlight the strengths and weaknesses of the models available in the literature, the authors built them and, after a validation process, performed a comparison. The analysis revealed the influence of the temperature at which the solid and fluid properties were computed as well as the need to consider the thermal losses and the solid thermal conductivity.
Starting from these outcomes, the authors developed a new and more complete model for the tank, namely the TES-PD model. The model accounts for the thermophysical properties' local time and space variability, and it embeds the thermal losses and the solid thermal conductivity. Features that avoid any uncertainty concerning the reference temperature definition, and considering the solid and fluid properties variability layerby-layer, allowed for a better performance prediction of both the storage and the plant in which it works.
In addition, the TES-PD model can take charge of real properties, paving the way for the study of innovative fluids, materials and configurations regarding thermal storage issues if these parameters are available from accurate measurements, dataset or CFD analysis. Finally, it is possible to claim that the performed analysis on the TES-PD model revealed that such differential and time-accurate description of the storage tank, even if one-dimensional, allowed a fast and effective prediction of the performance of both the tank and the storage plant; this is essential information for the preliminary design of innovative large-scale storage units operating with thermal storage.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Detailed Description of Most Used Numerical Models for Energy Storage Tanks
In this section, the most used numerical models for storage tanks are described with a high detail level. The logic of an increased complexity governs the description. Despite the models are already available in literature, only a few details were reported regarding their numerical implementation; thus, the authors present them in this section to highlight their peculiarities and the method that they developed for their mathematical resolution. The model m.1, presented for the first time in 1982 by Howell et al. [30], represents one of the most straightforward schemes for packed bed numerical description as demonstrated by its large use in the literature (see, e.g., [13,26,38,41,44,45]). The tank model consists of discretising the bed in N number of elements over the height of L; thus, ∆x = L/N represents the spacing within the bed layers. Inside any layers, the temperature and the pressure are assumed constant and uniform throughout the single bed element (the cross-sectional area is A). With this hypothesis, the pressure drop in the tank, ∆p, can be computed as: where C f is the friction coefficient, G =ṁ/A is the specific mass-flow rate,ṁ is the mass-flow rate, d is the particle equivalent diameter and ρ f is the fluid density. The fluid temperature distribution along the bed is computed as: The temperature of the bed, instead, holds as follows: Here, T f ,i is the temperature of the i-th fluid layer while T s,i is the temperature of the i-th bed's layer. The density of the solid material constituting the packed bed is denoted as ρ s , while c p,s and ε are the solid specific heat at constant pressure and the void fraction of the bed, respectively. The time-step, ∆t, is evaluated as ∆t = (ρ f L A ε)/(Nṁ). NTU is the Number of Transfer Units, U i is the overall heat loss coefficient, and ∆A i is the external area of a layer. The volumetric heat transfer coefficient, h v , depends on the thermal conductivity of the fluid, k f ; the hydraulic diameter, d; and the Nusselt number, Nu, according to the following expression: To compute such parameters, the correlations reported by Singh et al. [46] are employed; thus, the Nusselt number is evaluated as: Note that Re = (G d)/µ f is the Reynolds number, µ f is the dynamical fluid viscosity and Ψ = A s /A e is particle sphericity (i.e., the ratio between the surface area of a spherical particle of the same volume and the particle surface area). Finally, always according to Singh et al. [46], the friction coefficient is computed as: The above set of correlations is useful since they provide an algebraic description of storage tank and they are valid for packing material with a geometry different from the spherical one, only through the sphericity definition.
where and τ are, respectively, the lengths and the time scale defined as: refer to Desrues [42] for more details, while in the following, the adopted correlations are reported concerning spheres packing geometry.
• Microchannels structure. The friction factor is defined as: where α = L depth /L width defines the aspect ratio of the channels. The Nusselt number is defined as Nu = 7.541 (1 − 2.610α + 4.970α 2 − 5.119α 3 + 2.702α 4 − 0.548α 5 ) and it is used to compute the superficial heat transfer coefficient and the effective thermal conductivity: k s,e f f = (1 − ε) k s (A26) • Packing of spheres The hydraulic diameter, d h , is given as: while the hydraulic Reynolds, Re h , which is linked to the superficial one, can be defined as: and it is used to determine the friction factor based on the following equation: which represents the Ergun equation [49] for spheres packing, and the Wakao & Kaguei [50] correlation for Nusselt number Nu = 2 + 1.1(Pr) 1/3 (Re) 3/5 (A30) where Pr = (µ f c p, f )/k f is the Prandtl number and Re = (G d)/µ f is the Reynolds number.

Appendix B. Convergence Analysis and Grid Sensitivity
The section aims at performing a grid sensitivity analysis applied to the TES-PD model. In particular, a parametrical study is performed to tune the number of points required for the grid independence of the results. The analysis outcomes are depicted in Figure A1. In particular, Figure A1 shows the temperature trends inside the tank at the end of the first charge and discharge cycles, adopting an increasing number of subdivisions, N. Table A1 lists the main performance parameters. The reader can notice how the computational Time of Completion increases rapidly with a finer discretisation. The authors convey in tolerating a 2% variation of all the performance parameters, counting that further refinement of the grid leads to almost doubling calculation times, which affects the final use and the scope of a one-dimensional model. For these reasons, a number of layers equal to 200 are found a good compromise between numerical accuracy and computational efficiency. The same number of layers is used in the model m.3 basing on the formally equivalent numerical structure compared to TES-PD. x/L T/T max N=20 N=50 N=100 N=200 N=300 Figure A1. TES-PD model behaviour with increasingly refined grid level.