Development, validation and demonstration of a new Modelica pit thermal energy storage model for system simulation and optimization

Space heating applications account for a high share of global greenhouse gas emissions. To increase the renewable share of heat generation, seasonal thermal energy storage (STES) can be used to make thermal energy from fluctuating renewable sources available in times of high demand. A popular STES technology is pit thermal energy storage (PTES), where heat is stored underground, using water as a storage medium. To evaluate the use of PTES in an energy system, easily adaptable, publicly accessible and tool independent models are needed. In this paper, we improve an existing PTES model developed in the Modelica modeling language. The model is cross-compared with a more detailed and previously validated COMSOL model, considering different amounts of insulation, showing a deviation of 2–13% in the observed annual charged and discharged amount of heat. The results indicate that the presented model is well suited for early design stage and an exemplary case study is performed to demonstrate its applicability in a system context. Dimensions of system components are optimized for the levelized cost of heat (LCOH), both with and without subsidies, highlighting the importance of subsidies for the transition towards climate friendly heating solutions, as the gas boiler use is reduced from 47.6% to 2.7%.


Pit thermal energy storage systems for solar district heating
A large share of around 50% of the total energy demand in Europe is used for heating and cooling purposes (HRE 2019).As more than three-quarters of this demand is met by non-renewable energy sources, this sector is a large contributor to the production of greenhouse gas emissions (Eurostat 2022).Therefore, there is a large incentive in the European Union to increase the share of renewable heat sources in the heat supply scheme.A promising option in this regard is solar thermal energy from central heating plants or the utilization of unavoidable waste heat (Pelda et al. 2020).One of the main drawbacks of these technologies is the fluctuating energy supply over the year on a seasonal, weekly and daily basis (Gao et al. 2015;Hirvonen and Kosonen 2020).This is especially the case for the availability of solar heat during the summer and high heating demands during the winter (Bauer et al. 2010;Schmidt et al. 2004).This mismatch can be addressed by using thermal energy storages (TES).TES allow for the storage of thermal energy from sources which are not dispatchable on demand, such as solar or waste heat, and make it available at later times, bridging the gap between supply and demand (Dincer and Rosen 2011).This can maximize both the flexibility and performance of district heating (DH) systems (Guelpa and Verda 2019;Chicco et al. 2022) and enhances the capability to integrate renewable energy sources into DH systems (Sifnaios et al. 2023).However, there is a broad range of TES systems, and the selection of an appropriate system depends on many factors like the general application, operating conditions, storage period and economic factors (Ibrahim et al. 2008;Fournier et al. 2024).
Regarding their storage time, TES systems can generally be classified into long-and short-term storage systems (Dahash et al. 2021a).The concept of short-term storage implies charge and discharge periods of a few days and is often used, when the daily heat demand and the heat production vary.Long-term storage or seasonal thermal energy storage (STES) periods can instead range up to multiple months and are especially suitable for heat sources with seasonal variation in availability (Guelpa and Verda 2019).
STES that are either fully or partially covered by soil are also referred to as underground thermal energy storages (UTES).The most common types are: tank TES (TTES), pit TES (PTES), aquifer TES (ATES), mine TES (MTES) and borehole TES (BTES) (Heatstore 2021).This work focuses on PTES, which are characterized by high charging and discharging energy rates compared to ATES and BTES systems, which makes them especially useful for solar based DH (Dahash et al. 2020).A common limitation of central solar plants and PTES, however, is their dependence on the availability of large areas of free land at economically viable costs, limiting them so far mostly to DH grids in rural and suburban communities (Trier et al. 2018).
A PTES consists of a pit in the ground with inclined walls, that is enclosed with watertight liners (Fig. 1).One key point for the design of PTES is providing an efficient insulation to reduce heat losses.Typically, insulation is implemented by layers of glass wool, polyurethane and expanded polystyrene, which are attached to the watertight liners of the lid (Yang et al. 2021).A detailed technical description of PTES designs and insulation materials is given by Xiang et al. (Xiang et al. 2022).Factors like the operating temperature levels or the quality of thermal stratification of the storage fluid influence the performance of a PTES system considerably (Sifnaios et al. 2022).A good thermal stratification of the stored water can be achieved by measures like a proper diffuser design (i.e., inlet/outlet device), to reduce turbulences within the water and good thermal insulation at the top, to avoid cooling down of the uppermost water layer.Operating conditions are primarily defined by the integration in the overall heating network (Sifnaios et al. 2022;Pauschinger et al. 2018).This includes an aligned dimensioning of further system components and used materials, but also aspects like the hydraulic layout and control strategies (Pauschinger et al. 2018;Hesaraki et al. 2015).Because of the massive dimensions of PTES systems and the resulting demand for construction materials, all these aspects can have a high impact on the profitability of a project, resulting in a complex and cumbersome planning process.Thus, planners of such thermal systems tend to exploit numerical simulation tools that allow them to examine a wide range of boundary conditions producing the optimal TES for the considered site.Such planning tools assist in alleviating the risks emerging from the construction of such technologies and reduce investment cost.However, these planning tools require models with an accurate representation of PTES performance, which consider both the intricate design aspects and material choices as well as interdependencies with the overall system.Furthermore, flexibility and expendability are desirable features, as PTES are still a rather novel technology and therefore subject to continuous design changes.
To address these issues and support the further roll-out of PTES, this work presents a new open-source PTES model which is implemented in the modeling language Modelica (Olsson 2017).The model incorporates the most important design features of realized PTES systems and can be adapted to further developments of this technology easily, due to the modular and open approach of Modelica.A cross-validation study is carried out to assess the accuracy of the models in detail and finally a solar district heating (SDH) system with PTES is optimized in an exemplary study, to highlight a typical application and capabilities of the new model.

State-of-the-art modeling approaches for pit thermal energy storage systems for system simulation
For numerical simulations, a wide variety of tools exist that can be used to represent PTES performance and investigate design aspects.Such tools can be classified after their spatial discretization method into finite difference (FD), finite volume (FV) and finite element (FE).Alternatively, the tools can also be classified according to the level of detail involved in the simulation model and accordingly: component-level, technology-level and system-level.This choice plays a critical role in the required simulation time and thus the effort needed to accomplish the planning, especially the comparison of different designs and optimization, of such technologies.It is worthwhile to mention that the majority of TES models which are applicable for system simulations, are discretized in FD fashion.Thus, the following discussion pays attention to these aspects and addresses the recent advances.
To investigate integration of PTES on a systems level, with a focus on TES planning aspects, several simulation tools, such as TRNSYS (Solar Energy Laboratory 2021), Dymola (Systèmes 2020), or SimulationX (Iti 2021), are available.These allow us to gain insights into the interaction between the TES layout and the system.Many numerical models are developed as 1-D, 2-D or 3-D representations in TRNSYS (Xiang et al. 2022).For instance, the XST model (Type 342) and the ICEPIT model (Type 343) are typical models used for simulation of large-scale TES with cylinder or cone geometries, respectively.
Herein, Bai et al. (2020) developed a simplified underground STES that allows to simulate a TES with 1-D representation embedded in a 2-D soil model.The model is discretized using the FD method.It was calibrated against measured data from a PTES with a volume of 3000 m 3 , indicating acceptable agreement.Furthermore, Xie et al. (2021) improved the features of Type 343 (ICEPIT) and conducted a calibration study with the aid of measured data from the Dronninglund PTES with a volume of 60,000 m 3 .The results showed deviations between the measured and simulated annual amount of charged and discharged heat of below 5% and of below 10% for heat losses.
Alternative tools to TRNSYS are, e.g., Dymola and SimulationX, which are based on the Modelica modeling language.Using Modelica, several libraries have been developed to simulate energy systems.Such libraries can be free [Buildings (Wetter et al. 2014), AixLib (Müller et al. 2016), etc.] or commercial [e.g., TIL (Junior et al. 2008), ClaRa/ ClaRa + (Vojacek et al. 2023)].Compared to other modeling approaches, Modelica models show better reproducibility and expandability features that allow the users to deploy the codes in several Modelica-based simulation tools (e.g., SimulationX, OpenModelica, etc.).
The Buildings library (Wetter et al. 2014) is one of the most used libraries in both academia and industry for district energy simulations.It includes several tank models for simulating 1-D stratified storages.Yet, these models are restricted to residential applications (i.e., building energy systems).Moreover, they cannot represent all desired geometries (e.g., cones) and lack a proper soil model to simulate the heat transfer between the TES and the surrounding soil.
Thus, Dahash et al. (2022) further developed the stratified hot water tank model in the Buildings library.The work paid notable attention to extending the model capabilities for simulating UTES.It also developed a proper 2-D soil model that allows to host the TES in the soil.It is worth highlighting that both models are discretized using the FD method.The developed model was calibrated against a previously validated large-scale TES model for buried tanks and the outcomes underlined a good accuracy for simulating large-scale buried tanks, with a deviation in thermal losses between both model approaches of about 2%.Yet, the model is limited to simulating a cylindrical geometry.Thus, further efforts are required to extend its capabilities to represent other geometries (e.g., cones).
Most recently, Reisenbichler et al. (2022) further developed a model for large TES based on the hot water tank model embedded in the Buildings library.The developed model is further connected to a 2-D soil model discretized in an FD fashion.The model was calibrated using measured data from the Dronninglund PTES for the year 2015.The outcomes showed a discrepancy between the results and the measured data with a maximum deviation of 12% in the thermal losses.This deviation results from the fact that the model assumes the PTES to have a cylindrical shape, whereas the Dronninglund PTES is shaped like a truncated pyramid.Thus, several adjustments in the heat transfer coefficients for the lid, sidewalls and bottom were made to counterbalance the impact of the pit surface area.
Furthermore, Dahash et al. (2020) developed a versatile model that can simulate several TES shapes (e.g., cylinder, cone, pyramid, cuboid, or a combination) hosted in a 2-D/3-D soil model.The model is discretized in an FE fashion and simulated using COMSOL Multiphysics (COMSOL 2014).As validation is a cornerstone in the development of numerical models, the model was validated by comparing it against measured data from the Dronninglund PTES for the year 2015 (Dahash et al. 2021a).The validation outcomes highlighted a remarkable agreement between the measured data and the simulation results for both the 2D and 3D models with deviations in the charged and discharged amount of heat below 0.5%.The model can be used for component-level or technology-level simulations, but is limited in terms of system simulations due to several restrictions (e.g., simulation time, coupling to system model).
The foregoing discussion highlights the necessity to develop a PTES model that considers key design aspects, like pit geometry, accurately representing the thermal behavior of the storage volume and the surrounding ground, which can be incorporated in a system simulation environment.Furthermore, the model should be easily adaptable to new design aspects, publicly accessible and tool independent to promote a wide use.In this work, the authors therefore extend an existing Modelica PTES model and demonstrate its applicability for system design optimization within an exemplary study.

Numerical model development
The PTES model is defined with strict conformity to the Modelica standard, which is an open-source modeling language, actively developed by the Modelica Association (2024).It can therefore be used for simulation in a wide range of tools which support the Modelica language or even without such a software, by model export (Junghanns et al. 2021).Modelica is based on a concept called physical modeling, which denotes the acausal definition of physical laws by their equations, to promote reusability of components (Mattsson et al. 1998).Model components are typically collected and shared in libraries, which can range from low-level functionalities, such as thermal resistances or pipes (Tiller 2000;Zimmer et al. 2021), to sophisticated applications, like whole energy systems (Müller et al. 2016).Models of different libraries can be combined, allowing for the construction of high-level components by components from a low-level library.This is demonstrated by the reuse of simple model components for thermal conduction and convective processes from the Modelica standard library (Association 2020) to build up the new PTES model.
The original PTES model was described by Kirschstein and Formhals (Kirschstein 2020;Formhals 2022a).Since then, the model has undergone further development and has been released open-source in the Modelica Solar District Heating library [MoSDH (Formhals 2022b)].Therefore, only changes regarding the original model, which are already implemented in the MoSDH library, are presented in the following.
The pit volume model, which formerly divided the fluid volume in several layers of equal height, now divides the fluid volume in layers of equal volume.Since specific heat capacity and density of the storage medium are considered as constants only, this corresponds to layers of equal heat capacity.Between 30 and 90 °C the specific heat capacity and density of water vary by 0.5% and 3% (VDI 2010).
To calculate convective heat transfer due to thermal inversion, two new approaches known from the widely used Modelica libraries Buildings (Wetter et al. 2014) and Build-ingSystems (Nytsch-Geusen et al. 2013) have been implemented in addition to the originally implemented approach of the AixLib library (Müller et al. 2016).Each of these approaches has specific advantages and drawbacks, which have been discussed in detail in the original literature and previous studies (Zofer 2019).In both newly implemented approaches, an additional heat flow from bottom to top layers is used to compensate thermal inversion.For the Buildings approach, mixing effectiveness can be considered by means of a time constant τ, representing the mixing time (Eq.( 1)).In contrast to that, the BuildingSystems approach uses empirically tuned parameters G and n exp n to calcu- late the heat flow corresponding to the buoyancy induced mass flow rate using Eqs.( 2) and ( 3), where n denotes the number of volume elements, V the storage volume, ρ and c p medium density and specific heat capacity and ΔT the temperature difference of the adjoining volume elements: Another improvement concerns the modeling of the sloped walls of the PTES.While the heat transfer coefficient between storage medium and ground is calculated for the inclined slab, up to now, the inclination was not properly considered in the adjacent ground elements.Originally, elements from the step-shaped ground mesh were directly connected to the inclined pit elements and only heat transfer in horizontal direction was accounted for, whereas the vertical component was neglected.For the new approach, smaller elements are added adjacent to the pit wall, to better approximate the slope of the pit (see Fig. 2).Their cross-sectional area is defined equivalent to the area of the (1) corresponding elements of an inclined wall (triangular cross section).The heat capacity of such a ground element is located at its center of mass.The thermal resistance (R in Fig. 2) in the four directions (up, down, left, right) represents the respective effective resistance to the mass center in r or z coordinates.The upper and inner (left) heat ports are connected to the respective storage volume element, such that radial and vertical heat transport is now considered.

Validation study design
Validation represents the cornerstone in the development course of new dynamic models and tools.This is due to its importance in gaining trust in the developed model and it later allows the applicability of the model for planning and design.Consequently, this work includes a validation for the developed PTES model.This could be done using measured data from a real-world application for several years of operation.However, it is challenging to obtain such data of the required quality, for a technology as young as PTES.Furthermore, real-world UTES systems are difficult to implement accurately, since many parameters like soil properties are hard to determine with an accuracy higher than typical margins of errors for validation studies of about 5-10%, do not show homogeneous distribution of properties within the ground or change their values with time (e.g., saturation).Thus, the authors decided to cross-compare the developed Modelica model against a highly detailed and validated model, which was calibrated by fitting the thermal conductivity parameters of insulation material and surrounding ground to match the monitored behavior.
The model is described by Dahash et al. (2020), is implemented in COMSOL Multiphysics, using the FE method and can represent a wide range of TES geometries (cf.chapter 1.2).Readers are referred to original literature for an in-depth description of the model setup.The model is validated against measured data from a pit thermal energy storage with a volume of 60,000 m 3 , which is Dronninglund pit TES (DK).After parameter fitting the work showed a remarkable agreement between the simulated results and the measurements for the year 2015, resulting in deviations for the charged and discharged amount of heat as well as the amount of thermal losses of below 0.5%.Thus, this work uses the COMSOL model as a benchmark for the cross-comparison.
This work uses the cross-comparison framework described by Ochs et al. (2022) as basis, which was specifically designed for the comparison of different UTES models within the research project "giga_TES".The framework limits the comparison boundaries to the storage and, consequently, there is no need to include the entire DH system in the comparison.Yet, the characteristics of the DH system (e.g., heat supply, heat demand, supply temperature, return temperature, flow rate) play a dominant role in the charging and discharging of the storage integrated.Thus, the work developed a profile for the flowrate injected into TES to emulate the charging and discharging modes as shown in Fig. 3.The framework predefines the temperatures of 90 °C and 60 °C as DH supply and return temperature, respectively, as presented in Fig. 4. In this context, it is worthwhile to highlight that both foregoing figures demonstrate a seasonal operation for the storage comparison.
Furthermore, the framework (Dahash et al. 2020) describes a fully buried TES as the central comparison element with a volume of 100,000 m 3 , a height of 25 m and a slope of 30°.The TES has two inlet/outlet ports; one is installed at the top with 0.5 m underneath the cover and the other at 0.5 m above the bottom.As the TES is fully hosted in the soil, the framework defines a depth of 30 m for the soil below the TES bottom.Besides, the heat transfer coefficient between the ground surface and ambient air is set to 25 W/ (m 2 K).In this regard, Fig. 5 shows the TES hosted in the soil domain with the most relevant boundary conditions and dimensions.The surrounding soil has the following thermal properties: density of 2000 kg/m 3 , specific heat capacity of 880 J/(kg k) and thermal conductivity of 1.5 W/(m K).The initial soil temperature is set at 10 °C, which is equivalent to the initial TES temperature.
The ambient temperature is expressed as below: Due to the lid interaction with ambient air and the elevated temperature at TES top, it is important to compare the TES under several lid insulation cases.Even though (4) experience so far showed that such large-scale storages can be built with no insulation on the sidewalls and bottom to reduce capital costs as the cost-to-benefit ratio showed low gain, recent studies proposed that insulation can be beneficial under certain boundary conditions (Dahash et al. 2021b).With increasing heat prices higher investment costs for insulation might become feasible, depending on economical boundary conditions.Thus, the framework divides the comparison into two primary cases: "non-insulated" and "insulated".Each item is further divided into three minor cases considering different insulation quality of the TES lid.In this regard, the noninsulated cases are labeled with (Case 1), whereas the insulated cases are assigned to (Case 2).Table 1 reports a summary of the cases considered in the comparison and their corresponding characteristics.
The simulations are run for 10 years for the uninsulated TES cases (case 1), while cases 2 (insulated) are run for 5 years.This is attributed to the fact that soil has an initial temperature of 10 °C and gradually heats up during initial storage cycles until an equilibrium is reached on a seasonal scale.After this initial phase, operation conditions between succeeding storage cycles are virtually identical, resulting in no further gain of information.This preheating phase might last between 3 and 5 years depending on the insulation quality, TES size and geometry and operation temperature range.For the insulated cases (case 2), less heat is injected into the surrounding ground resulting in a smaller thermally affected zone and a quicker equilibrium.Thus, the underlying validation framework defined simulation times of 10 years for case 1 and only 5 years for case 2 and pinpointed that the results analysis should be performed for the last year (10th year for case 1 and 5th year for case 2).
The results analysis should be performed for outcomes in an hourly resolution.In case further info is needed on the cross-comparison framework, the reader is referred to Dahash et al. (2020).

Exemplary application study design
To demonstrate the functionality of the developed PTES model for system simulations, an SDH design optimization study is carried out as an exemplary application, using a fictional scenario mostly based on German boundary conditions.The main system component dimensions are subject to optimization, to result in the lowest possible levelized cost of heat (LCOH), both with and without subsidies.

System model
Figure 6 shows a block diagram sketch of the system model, which is based on models from the MoSDH library (Formhals 2022a).It consists of a heat load connected to a DH grid that is supplied by a solar thermal collector field, a PTES, a heat pump and a gas boiler.Both heat load and grid are aggregated into a single component each.In this example, the size of the solar collector area per heat demand, the pit storage volume per collector area and the heat pump capacity per storage volume are subject to optimization.
The heat demand of the used heat load curve amounts to a total of 25 GWh/a and is constructed by adaption of monitoring data from a university DH grid to test reference year conditions in Darmstadt, Germany (Formhals et al. 2021).The heat load is met by cooling the hot heat carrier from the grid supply line to a return temperature of 30 °C.The grid is represented by a single component to emulate buried pre-insulated pipes (DN250, reinforced insulation) with a depth of 1.2 m and a length of 400 m divided into 10 segments.Soil parameters were chosen according to the benchmark test case scenario (cf.Sect."Validation study design").
A return flow mixing component is placed between the DH grid and the heat generation and storage components, which is used to lower down the supply temperature, by mixing with volume flow from the return line, if source temperature exceeds the temperature setpoint.This setpoint is determined according to a heating curve that is reliant on the ambient temperature and varies between 70 °C and 80 °C.

Fig. 6 Simulation model for the application example
The solar thermal collector field component models an array of collectors consisting of 5 collectors in series with an inclination angle of 35°, while the number of rows is varied by the optimization scheme.The modules are parametrized according to manufacturer data for a collector for ground-mounted flat plate collectors (SPF-Institut für Solartechnik 2018).Since the solar collectors use a 20% glycol mixture as heating medium, they are separated to the DH system by a heat exchanger with an exchange surface area of 0.2 m 2 per square meter collector area and a heat exchange value of 1500 W/(m K).Following the instantaneous heat load, the heat supplied by the solar field is either fed into the DH grid, into the PTES or a combination of both.The collector field is switched off if the PTES temperature exceeds 90 °C.
The PTES has a fixed depth of 10 m and the walls have a slope of 0.5, while the volume is subject to optimization.The storage pit is modeled by 15 layers of equal volume.Ground parameters are chosen according to the benchmark test case (cf.Chapter 2.2).No insulation is used for walls and bottom and the lid consists of 3 layers of polyethylene of 10 cm thickness and a thermal conductivity of 0.03 W/(m K) each.The mid-diffusor is placed at the volume center at a height of 6.8 m.
The heat pump model is parametrized according to manufacturer data of a hightemperature heat pump.The maximum heating power map is linearly scaled to meet the heating power defined by the optimization algorithm in the operation point W30/ W80, while the COP map is kept unchanged.The heat pump's source side is connected to the mid port of the PTES and it can feed into the DH grid and the top of the storage.The load supply temperature is set to 3 K above the reference supply temperature of the DH grid.During the heating season between October and April the heat pump is switched on, if the storage temperature 1.5 m below the lid falls 5 K below the reference supply temperature and the bottom diffusor temperature is above 23 °C.It is turned off, when the respective storage temperature exceeds the reference supply temperature, the storage bottom temperature falls below 18 °C or during summer season.
The gas boiler has a maximum thermal power of 10 MW and can therefore supply all of the required heat demand and, accordingly, serves as a back-up heating unit.However, it is only turned on, if the remaining heat sources cannot supply all the heat demand on the required temperature level.When the temperature at the top of the PTES falls 3 K below the reference supply temperature, the boiler is used to shift the temperature up to the desired level by mixing.Only when the storage temperature falls 10 K below the setpoint, no heat is extracted from the storage.The efficiency of the gas boiler is defined according to an existing boiler at the TU Darmstadt university campus of comparable size and varies between 89.5% and 93%.

Economic and environmental optimization
LCOH is calculated after Eq. ( 5), with maintenance costs M a , fuel costs F a , heat delivery Q a and return rate r.The different lifetime of system components is considered using annuities AN (Eq.6), evenly distributing investments costs over the lifetime a life of each component: Table 2 summarizes the used investment cost functions, annual maintenance costs and expected lifetimes, which are taken from literature and adapted by inflation for 2021 (Destatis 2022).
Table 3 shows the used fuel costs and subsidies.Electricity costs correspond to the annual average for industrial electricity in 2021.Since gas prices are still in a recovery phase from the price-shock due to the war in Ukraine and future projections strongly differ based on scenario and location (Schlund et al. 2023), a notional gas price of 10 € − ct/kWh was assumed in this study.Currently, many countries subsidize the construction and operation of renewable energies for DH.Exemplary, the new German funding scheme for efficient DH systems "BEW" [Bundesförderung Effiziente Wärmenetze (BMWK 2022)] was applied, in which up to 40% of the investment costs can be funded, up to a maximum of 100 Mio.€.Additionally, the funding is limited to the profitability gap, which is required to bring the business case to break-even.Since calculation of this gap depends on the specific case, it was not considered for this generic study due to simplification.Furthermore, heat from solar thermal collectors is supported by 1 € − ct/kWh and heat pump operation is subsidized to a maximum of 9.2 € − ct/kWh depending on the seasonal performance factor (SPF).
Optimization is carried out using the NMinimize function of Mathematica (Wolram Research 2021).The Nelder-Mead method is chosen (Nelder and Mead 1965), a direct search method for nonlinear optimization, often called simplex-method.
(5)  Heat pump operation subsidy min((5.5 − (6.8-17/SPF)*0.75)*(SPF/(SPF-1)),9.2) BMWK ( 2022) Equation ( 7) shows the target function of the optimization, which is the LCOH after Eq. ( 5).Optimization variables are defined as relative numbers, to account for the interdependence of solar aperture area A sol , storage volume V pit , nominal heat pump capacity P th,HP and annual heat demand Q a .Definition of the key design parameters as relative sizes allows for better comparability of systems with different total heat demand (Mauthner and Herkel 2016).Constraints for optimization variables are defined as relatives accordingly (Eqs.8-10): Calculation of the LCOH is carried out in the system model (LCA component in Fig. 6) during runtime of the simulation, according to the equations given in this chapter.The overall optimization algorithm is run in Mathematica, which is coupled to the simulation environment SimulationX by the COM interface (Iti 2021; Wolram Research 2021) to setup the simulation models, start simulations and receive the results of each iteration.

Model cross-validation
In order to justify the applicability of the model for planning purposes, it is crucial to first verify its results.Thus, the model was cross-compared against a validated model with the aid of the comparison framework described in Sect."Validation study design".The following subsections report the results of the comparison framework.For the sake of simplicity, the focus is set on two selected cases, whereby one case is a PTES without insulation on the side walls and bottom (case 1b) and the other has insulation (case 2b).Both cases have the same insulation on the top lid (Table 1).The other 4 cases, which differ only in the amount of insulation on the top, are summarized and will be briefly reported.

Energy balance
Herein, Fig. 7 reveals the total thermal losses and its breakdown for two selected cases: case 1b which assumes no thermal insulation on the pit walls and case 2b, which considers insulation to the ground.For case 1b, the simulations ran for 10 years allowing the PTES to operate in a quasi-steady-state post to passing the preheating phase in which the ground temperature drastically increases due to the difference between the (7) min TES operation temperature and ground temperature.In contrast to that, case 2b is only simulated over five years, since less heat is put into the ground and consequently, steady operation is achieved earlier (cf.chapter 2.2).Obviously, Fig. 7 emphasizes a deviation in the total thermal losses for both selected cases between the two compared models.For case 1b, the COMSOL TES model reports total thermal losses of 1208 MWh, whereas the Modelica PTES model has total thermal losses of 1322 MWh producing a deviation of ~ 11%.It is observed that the side losses contribution accounts for most of the aforementioned deviation between both models.On a brighter note, the deviation notably decreases down to 0.5% for the total thermal losses in case 2b.Yet, a deviation for the side contribution can still be observed, which is counterbalanced by the top thermal losses in this case.
Due to the irreversible thermal losses experienced by TES systems, the amount of the extracted thermal energy is less than that injected in the charging phase.Thus, it is also of importance to consider the charged and discharged energy in the comparison framework.Subsequently, Fig. 8 presents a comparison of charging/discharging energy for both selected cases.Therein, the positive values represent charged energy and negative values for discharged energy from PTES.For case 1b, the deviation in the charged energy amounts up to 7% and increases to 13% for the discharged energy.For case 2b, the observed deviations for both charging and discharging fall below 7%, indicating a better match between the Modelica and the COMSOL model.
Table 4 summarizes the main results of the remaining 4 cases.The main conclusions from the previous paragraphs detained by case 1b and 2b still hold true: -Overall, heat losses are smaller for the insulated cases and a smaller deviation to the COMSOL model is observed as well.Table 4 Energy losses through storage bottom (Q bot ), side (Q side ), top (Q top ) and overall (Q loss ) as well as charged and discharged amount of heat (Q ch/dis ) for the COMSOL and Modelica models (excluding cases 1b and 2b, which are described in detail above) -The highest relative deviation to the COMSOL model is observed for thermal losses through the storage bottom, but since the absolute amount is small compared to the losses through the wall and the top, this deviation is negligible for the overall performance.-The amount of charged and discharged energy is generally lower in the Modelica models with a range of − 2% to − 13%.

PTES temperature
Figure 9 compares the water stratification profile inside the PTES at selected relative heights of 0.05, 0.1, 0.25, 0.5, 0.75, 0.9 and 0.95 for case 1b, where a relative height of 0.95 represents a probe at 95% height of the storage volume, i.e., close to the PTES lid.It can Fig. 9 Comparison of temperature profile at selected relative heights (h*) of the TES for uninsulated case 1b at the 10th year be seen that the fluid temperatures in the top of the storage of both models diverge right after changes between charging to discharging mode, but converge to similar values for longer operation under steady conditions.For fluid temperatures closer to the bottom however, the Modelica model shows an underestimation of the maximum observed temperatures during charging, with a good agreement for temperatures during discharging.
Figure 10 shows the temperature distribution of the storage fluid of both model options for scenario 2b.The overall characteristics of the curves and deviations are the same as in case 1b.However, overall, higher temperatures can be observed in this case in comparison to case 1b.

Exemplary study results
The exemplary study described in Sect."System model" serves to highlight a typical application for the PTES model, which is the determination of central component

Component dimensions of the optimal systems
The optimization algorithm converges after 49 iterations for the scenario without subsidies and 46 iterations for the scenario with subsidies (Fig. 11).Consideration of subsides results in a significantly larger size of the solar thermal collector field with a relative aperture area per annual heat demand of 2.75 m 2 /MWh, whereas the optimal system without subsides results in a size of 1.3 m 2 /MWh.In contrast to that, the relative PTES volume per aperture area is comparable for both scenarios, the system without subsides having a specific volume of 3.17 m 3 /m 2 and the system with subsides of 3.07 m 3 /m 2 .In the case of the heat pump, the consideration of subsidies has a comparable effect as for the solar collector area, with a heating power approximately twice as high for the subsidized system (0.04 kW/m 3 and 0.08 kW/m 3 , respectively).

Shares of the overall energy supply for the optimized systems
Figure 12 shows the overall share of the heat supply of the system components for each simulated model.It is obvious, that consideration of subsides results in a much higher share of solar thermal energy, PTES as well as heat pump and accordingly a lower share of the gas boiler.The latter has a share of 47.6% without subsidies and only a remaining share of 2.7% with subsidies.Overall, only a small share of the solar thermal yield is utilized directly without storage, amounting to 11.7% for the case without subsides Fig. 11 Evolution of the optimization variables A sol /Q a (aperture area per annual heat demand) V PTES /A sol (storage volume per aperture area) and P th.HP /V pit (heat pump power per storage volume) over the optimization algorithm iterations with and without subsidies Fig. 12 Evolution of the share of energy supplied by the system components over the optimization algorithm iterations with and without subsidies and 8.4% for the case with subsidies.Unsurprisingly, the amount of thermal energy supplied by the PTES and the heat pump are much larger if subsidies are applied (73.1% and 16.1%) in contrast to the scenario without subsidies (38.1% and 2.9%).

Discussion
The cross-comparison of the Modelica and the COMSOL storage models reveals significant differences in the calculated losses to the ground, with generally higher deviations for the models which do not consider insulation.This indicates a better representation of thermal conduction processes with larger temperature gradients by the COMSOL model, which is not surprising, considering the finer and more sophisticated meshing scheme of the FEM software.This observation is supported by the fact, that the deviation between both models' relative losses is highest at the storage bottom, where the meshing scheme of the Modelica models results in the largest element sizes.However, since losses from PTES systems are generally highest at the top, where high storage temperatures prevail, absolute deviations are higher at these areas as well.Consequently, the lower accuracy of the Modelica models in calculating bottom losses is negligible for the overall deviation.Nevertheless, the overall values for storage losses as well as charged and discharged heat amounts still show a deviation by the Modelica model of around 5-10%.While this deviation indicates potential for improvement of the model, it is still in a range of other uncertainties like cost estimates, which are generally deemed acceptable during the concept stage of a project.It can be well explained by the fundamentally different levels of detail and associated numerical efforts of both approaches.The COM-SOL model has a much finer discretization, both within the storage volume as well as the ground.As a result, the Modelica models overestimate the amount of thermal mixing within the storage.This can be observed by the attenuated storage layer temperature profiles.The small number of modeled layers within the Modelica model is not able to fully account for a sharp stratification profile, which is a small mixing zone between the hot and the cold fluid within the storage.Furthermore, the simple FD meshing scheme of the Modelica model for the ground, does not allow for a good representation of the inclined boundary between storage volume and ground.This is an inherent drawback of the FD scheme, which was partly tackled by the introduction of smaller elements at the boundary (cf.chapter 2.1).Future work should focus on this aspect to improve the model accuracy.
Moreover, the accuracy of the Modelica model has to be assessed with consideration of its practical field of application.While the COMSOL model is a more accurate representation of the storage itself, it is not suited for extensive studies coupled to a dynamic energy system model.The Modelica model, however, is dedicated for system simulation studies, mainly for pre-design applications.The main purpose of these pre-design studies is to determine system component dimensions, which are likely to result in the best possible system performance later on.Input parameters like prices, energy demand, solar yield or thermophysical ground parameters have a high uncertainty at this planning stage.Consequently, the accuracy of the PTES model should be compared to these deviations between planning and actual values as well.For this reason, it is less important to improve the accuracy of the Modelica model by increasing its level of discretization, than being able to do a large number of simulation runs and include a wide range of parameter values.For the given reasons, the accuracy of the model can be regarded as suitable for system simulation.For a real-world application, the Modelica model would be used within a system model to determine the dimensions of the main system components during pre-design, whereas the COMSOL model would applied in a later stage, to optimize the storage design and assess thermal effects on the surrounding groundwater.
However, the presented results indicate that there is room for improvement of the Modelica model.Following the discussed inaccuracy for larger temperature gradients, both inside the storage volume as well as within the ground, the mesh discretization scheme should be improved.An obvious option would be the decoupling of the storage volume and ground discretization, to allow for a larger number of storage layers, without a significant impairment of the numerical effort.Furthermore, discretization within the ground should be improved within the boundary region to the storage volume, to reduce the discussed inability to account for large temperature gradients within the ground.It should be noted, that groundwater flow cannot be accounted for in the Modelica model and was not subject of this study.
The shown application example study primarily proves the applicability of the developed model for system pre-design.Optimization of main component dimensions with the aim of minimizing the LCOH is a typical problem during a projects pre-design stage.During this stage it is important to obtain a preliminary design, which will be refined during the following project stages.In such a later stage, a detailed storage model like the COMSOL model should be used to get an in-depth understanding of the storage characteristics.Without going into the results of the application example in detail, it is obvious, that subsidies are currently still required to result in a SDH system with both a high share of renewable energy as well as a competitive LCOH.This might change if gas prices increase significantly or CO 2 certificate costs are considered.Both would make gas-based heat less attractive and consequently lead to a lower fossil share.Nevertheless, under the assumed economic and subsidy boundary conditions and without setting a target for the renewable share, the optimization results in a system with only a small share of the auxiliary fossil boiler, by minimization of costs.Consequently, the funding scheme succeeds in setting the right boundary conditions and incentives.
A further key message of the application study is the importance of a system view for the design of storage components.Their behavior is completely defined by the interactions with heat generation and demand.Even though the system with subsidies has a solar collector area about twice the size of the system without subsidies, the amount of solar energy which is directly used, is even a little bit smaller.This is in line with a general observation for SDH systems, where solar thermal collector efficiency is higher for systems with lower solar fractions.

Conclusions
In this study, improvements on an existing PTES model in Modelica are reported, the model is cross-compared with a validated COMSOL model and applied in an energy system optimization study.The following conclusions can be drawn: • The presented model shows a deviation of 5-10% regarding storage losses and charged/discharged heat in all considered cases compared to the COMSOL model.
• The deviation is noticeably lower for the insulated cases.
• The model is well suited for energy system optimization in an early planning stage where it is simulated together with other system components to derive the component dimensions.• The COMSOL model is more detailed and produces more accurate results, but is more suited for later project phases, where the focus is set on individual component design and thermal impact on the groundwater, due to its significantly higher numerical effort and possibilities for coupling to other system components.• Under the assumed boundary conditions (e.g., gas price, no emission penalty) subsidies are still needed to reach a high share of renewables in district heating as a result of a purely cost-driven evaluation.

Fig. 1
Fig. 1 Typical design of a PTES with three diffusers (inlets/outlets)

Fig. 2
Fig. 2 Model mesh.One ground element connected to the storage volume is depicted in detail with thermal resistances R and heat capacities.The different soil colors refer to the ground and the embankment, respectively

Fig. 7
Fig. 7 Breakdown of thermal losses for case 1b and case 2b for the compared models

Fig. 8
Fig. 8 Charged and discharged energy for case 1b and case 2b for the compared models

Fig. 10
Fig. 10 Comparison of temperature profile at selected relative heights (h*) of the TES for insulated case 2b at the 5th year

Table 2
Invest I 0 and maintenance cost M a relations for the used components

Table 3
Fuel prices and subsidies used for the optimization study