DOI: 10.3303/CET23103145 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Paper Received: 12 May 2023; Revised: 05 June 2023; Accepted: 18 June 2023 
Please cite this article as: Zhang S., Ocłoń P., Yildirim M., Varbanov P.S., Arsenyeva O., Kapustenko P., 2023, Validation of the Coupled Heat 
and Moisture in the Soil for Underground Thermal Energy Storage Systems, Chemical Engineering Transactions, 103, 865-870  
DOI:10.3303/CET23103145 
  

 CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 103, 2023 

A publication of 

 

The Italian Association 
of Chemical Engineering 
Online at www.cetjournal.it 

Guest Editors: Petar S. Varbanov, Panos Seferlis, Yee Van Fan, Athanasios Papadopoulos 

Copyright © 2023, AIDIC Servizi S.r.l. 

ISBN 979-12-81206-02-1; ISSN 2283-9216 

Validation of the Coupled Heat and Moisture in the Soil for 

Underground Thermal Energy Storage Systems 

Sheng Zhanga,*, Paweł Ocłońa, Mehmet Yildirima, Petar Sabev Varbanovb, Olga 

Arsenyevab, Petro Kapustenkob 

aEnergy Department, Faculty of Environmental and Energy Engineering, Cracow University of Technology, Al. Jana Pawła II  

 37, 31-864, Cracow, Poland 
bSustainable Process Integration Laboratory – SPIL, NETME Centre, Faculty of Mechanical Engineering, Brno University of  

 Technology – VUT Brno, Technická 2896/2, 616 69 Brno, Czech Republic  

 sheng.zhang@doktorant.pk.edu.pl 

The Renewable Energy System for Residential Building Heating and Electricity Production (RESHeat ) system 

is a typical application of the utilisation of sun-tracked Photovoltaic Thermal (PVT) panels and underground heat 

storage units. This work introduces a MATLAB simulation model coupled with moisture and heat transfer 

processes. Validation of the proposed model from experimental measurements was conducted. The simulations 

were performed for a single year of RESHeat system operation for demo sites located in Cracow City. Results 

showed that the relative error between calculation and measurement varies from 0.258 % to 5.829 %, along 

with the same temperature trend. A Control Volume method is used to simulate coupled heat transfer in the 

ground. 

1. Introduction 

Renewable Energy (RE) makes great contributions to the natural environment and sustainable development. It 

provides a solution for solving the problem of excessive emission of Greenhouse Gases (GHG) by large 

traditional fossil fuels consumption. Solar energy resource, as a promising and reliable alternative renewable, 

is widely used in the architectural and industrial sectors. However, the biggest obstacle in applying solar energy 

as a RE lies in the fact that there is a mismatch between the highest availability and highest demand influenced 

by the natural climatic condition. The Renewable Energy System for Residential Building Heating and Electricity 

Production (RESHeat) passes a novel, fully renewable energy source on integrating heating, cooling and 

electricity production with the Underground Heat Storage Unit (UHSU). The basic principle of the proposed 

system is using solar energy as a primary RE. A series of sun-tracked PVT panels and sun-tracked Solar 

Collectors (SCs) are arranged to generate heat and electricity for thermal utilisation and power generation. A 

ground source heat pump serves the heating and cooling of the building. The UHSU exemplifies the novelty of 

gathering thermal energy from the waste heat and excess heat from SCs and cooling processes of PVT during 

the summertime and delivers heat for Domestic Hot Water (DHW) use by the hot water storage tank (an 

insulated tank) whilst regenerating the ground soil to avoid the decrease of heat pump efficiency via a heat 

accumulator (a non-insulated tank) (Ocłoń et al., 2023). 

Underground Heat Storage Units are widely recommended as an excellent approach for designing the energy 

framework in civic buildings. Ocłoń et al. (2020) investigated a model of the radiator used for PV panel cooling 

by the 3D control volume method. It is displayed that lowering the PV panel’s temperature could increase the 

gross efficiency of sun-tracked PVT units. Ocłoń et al. (2021) introduced an optimisation method for a heat pump 

system coupled with the UHSU to obtain the maximum share of renewable energy sources. The Jaya algorithm 

could lead to a total of 93 % heat energy which was covered by utilising UHSUs. Xu et al. (2022) presented a 

cylindrical model of seasonal thermal storage systems with solar collectors. Daneshazarian et al. (2022) 

developed a two-step optimisation model to match the thermal storage medium (TSM) in UHSU for buildings. 

Zhang et al. (2022) provided a comprehensive review of the application of hybrid renewable energy systems in 

865



buildings including PCM thermal energy storage. Li et al. (2023) developed a nearly net-carbon refrigeration 

approach with an underground cold storage unit. For the selected area, this promising system could fulfill the 

indoor cooling requirements and showed a great potential for future buildings own to the greener characteristic. 

This paper introduces a coupled moisture-heat model of an underground heat storage unit in RESHeat. The 2D 

model concerns solar-assisted ground resource heat pump, the PVT, solar collectors and UHSU which takes 

the influence of ground heat transfer into account. The simulations were performed for a single year of RESHeat 

system operation for demo sites located in Cracow City. A Control Volume method under the MATLAB platform 

is used to simulate coupled heat transfer in the ground. The coupled model is verified by the actual 

measurements of the demo site in Poland. 

2. Mathematical model of the coupled model 

The proposed heating system consists of sun-tracked PV panels and solar collectors, water type ground heat 

pump and an underground storage system (Figure 1). It utilises the thermal energy from PV panels cooling and 

SCs while waste heat and excess heat from PV panels and SCs is accumulated in the ground via the UHSU 

and returned for additional thermal energy preparation during heating seasons. 

 

Figure 1: Schematic description of the proposed design 

  

Figure 2: Geometric scheme for calculation of the ground temperature 

The geometric of the underground heat storage tank is followed in Figure 2. It assumes the symmetrical partition 

for simplifying the heat transfer for the cross-section of the calculating domain. The tank burial depth 𝐻1 is 

assumed as 5 m, and the width of the domain 𝑊 is assumed as 12 m. The total height of the domain 𝐻1 + 𝐻2 is 

20 m. The outer radius of the tank is equal to 1.5 m. The mathematical model of underground tanks is developed 

in MATLAB. The Finite Volume Formulation is used for defining a ground domain with all the relevant boundary 

conditions. In Cracow, the yearly average ground temperature is set to 7.5 °C. With the symmetrical condition, 

the left side of the domain is insulated. In addition, the right side of the calculation domain is insulated as well. 

2.1 Transport equations 

Moisture condition contains liquid water and vapors accompanied by soils underground. Liquid water flow is 

displayed by Richards’ equation, which represents the movement of water in unsaturated soils. In this equation, 

water flow 𝑞𝑤 is proportional to the product of hydraulic conductivity 𝐾 and gradient in water potential, where 

matric potential 𝜓 and local gravitational potential 𝑔𝑧 are considered. The relative location 𝑧 is measured from 

the datum plane and gravitational acceleration 𝑔 can be considered constant. 

𝑞𝑤 = −𝐾∇(𝜓 + 𝑔𝑧) (1) 

866



Vapor flow 𝑞𝑣  is derived based on Fick’s law which describes the relationship between diffusion flux and 

concentration. Here water vapor flow is driven by a gradient in vapor concentration, which depends on the 

temperature 𝑇 and humidity ℎ. The vapor flow can be calculated as the sum of temperature-driven diffusion flux 

𝑞𝑣,𝑇 and isothermal diffusion flux 𝑞𝑣,𝑖. 

𝑞𝑣 = 𝑞𝑣,𝑖 + 𝑞𝑣,𝑇 = −𝐷𝑣𝑐𝑣
′ ∇ℎ − 𝐷𝑣ℎ𝑠∇𝑇 (2) 

Where 𝐷𝑣 is vapor diffusivity, 𝑐𝑣
′  saturated vapor concentration, 𝑠 the slope of saturated vapor concentration 

function curve. 

Heat flow underground consists of sensible heat component, latent heat component and the convective thermal 

component brought by flowing water. Sensible heat flow is formulated as the Fourier’s law of heat conduction, 

where the thermal energy flow is calculated as the product of thermal conductivity 𝜆  and the gradient in 

temperature. Latent heat flow is related to the phase transition process of water vapor, accompanied by thermal 

energy absorbed and released into soil. It can be derived as the sum of latent heat of vaporization 𝐿 and heat 

energy in water 𝑇𝐶𝑤, where 𝐶𝑤 is the specific heat of water. 

𝑞ℎ = −𝜆∇𝑇 + (𝐿 + 𝑇𝐶𝑤)𝑞𝑣 + 𝑇𝐶𝑤𝑞𝑤 (3) 

The continuity equation for the mass in this closed underground system states the relationship of water flow 

including liquid water and vapor, and the change of water content per volume 𝜌𝑤𝜃, where 𝜌𝑤 is the density of 

water and 𝜃 the volumetric water content. 

𝜌𝑤

𝜕𝜃

𝜕𝑡
= −∇ ∙ (𝑞𝑤 + 𝑞𝑣) (4) 

The law of conservation of energy provides the relationship between the change of temperature and thermal 

energy flow. In the underground thermal system, specific heat is composed of three parts, the heat in liquid 

water, 𝐶𝑤𝜌𝑤𝜃, the heat stored in the soil, 𝐶𝑚𝜌𝑚𝜒𝑚, and the heat from the air fraction, 𝐶𝑎𝜌𝑎𝜒𝑎. 𝐶𝑚 and 𝐶𝑎 are the 

specific heat of soil component and air part respectively. 𝜌𝑚 represents the density of soil and 𝜌𝑎, the density of 

air. 𝜒𝑚 and 𝜒𝑎 are the volumetric fractions of solid minerals and air. The change of specific heat in air part makes 

little contribution to the energy flow compared to the changes in liquid and solid components, own to the amount 

of specific heat stored in the air is so small (Kroener et al., 2014). Therefore, it can be neglected when calculating 

the thermal energy. 

(𝐶𝑚𝜌𝑚𝜒𝑚 + 𝐶𝑤𝜌𝑤𝜃)
𝜕𝑇

𝜕𝑡
= −∇ ∙ 𝑞ℎ (5) 

2.2 Parameter functions 

Volumetric water content and hydraulic conductivity are significant properties for evaluating the underground 

moisture degree. Campbell G.S. displayed the relationship between hydraulic matric potential, volumetric water 

content, and hydraulic conductivity. 

𝜓𝑒 is the air entry potential, 𝜃𝑠 the saturated volumetric water content, 𝐾𝑠 the saturated hydraulic conductivity. 

And 𝑏 is the shape factor. The saturated volumetric water content is the maximum volume of water per unit 

volume of soil. 

𝜃 = {
𝜃𝑠(

𝜓𝑒

𝜓
)1/𝑏 ,   𝜓 < 𝜓𝑒

𝜃𝑠,                 𝜓 ≥ 𝜓𝑒

 (6) 

𝐾 = {
𝐾𝑠(

𝜓𝑒

𝜓
)2+3/𝑏 , 𝜓 < 𝜓𝑒

𝐾𝑠 ,                   𝜓 ≥ 𝜓𝑒

 (7) 

Humidity is the concentration of water vapor existed in the air. It is equal to the fraction of vapor concentration 

𝑐𝑣 and the saturated vapor concentration 𝑐𝑣
′ . In addition, it is related to matric potential and temperature. Here, 

𝑀𝑤 is the water molecular weight, 𝑅 the gas constant, and 𝑇𝑘 is temperature in Kelvin. 

ℎ = exp (
𝑀𝑤𝜓

𝑅𝑇𝑘
) (8) 

The saturation vapor pressure curve separates the two-phase state and the superheated vapor state in the 

temperature-entropy chart. The term 𝑀𝑤/(𝑅𝑇𝑘) functions the convertor from pressure to concentration in the 

867



saturation vapor concentration curve. Therefore, it can be calculated for the slope of saturation vapor 

concentration curve 𝑠, where Δ is the slope of the saturated vapor pressure curve. 

𝑠 =
∆𝑀𝑤

𝑅𝑇𝐾
 (9) 

In the vapor pressure curve, the Arden Buck equation provides a series of empirical correlations which combine 

the saturation vapor pressure and temperature of moisture air. The temperature is in Celsius. 

𝑒𝑠(𝑇) = 0.611exp (
17.27𝑇

𝑇 + 237.3
) (10) 

Here, the saturated vapor concentration is calculated as: 

𝑐𝑣
′ =

𝑒𝑠𝑀𝑤

𝑅𝑇𝐾
 (11) 

The rate of evaporation of water into the atmosphere can determine the mutual coefficient of diffusion of water 

vapor and air. Water vapor diffusivity 𝐷𝑣 is related to the air in the porosity and the binary diffusion coefficient of 

water vapor 𝐷0(𝑇, 𝑃), which depends on temperature and pressure. 

𝐷𝑣 = 𝐷0(𝑇, 𝑃)𝜀(𝜒𝑎) (12) 

𝜀(𝜒𝑎) = 𝛽(𝜒𝑎)𝑚 (13) 

The factor 𝛽 is 0.9 and 𝑚 is 2.3. The binary diffusion coefficient of water vapor in atmosphere is shown as: 

𝐷0(𝑇𝑘 , 𝑃) = 𝐷0(237.15K,101.3kPa) (
𝑇𝑘

273.15K
)
1.75

(
101.3kPa

𝑃
)      (14) 

For the standard measurement condition, 𝐷0(237.15K,101.3kPa) = 2.12 × 10-5 m2 ∙ s-1. 

2.3 Partial differential equations 

The transport equations in 2.1 and parameter functions in 2.2 above lay the function of partial differential 

equations for mass conservation and energy conservation. In Eq(15) and Eq(16), there are 2 unknown 

parameters to be solved, matric potential and temperature. 

𝜌𝑤

𝜕

𝜕𝑡
𝜃(𝜓) = 𝛻 ∙ (𝐾(𝜓)𝛻(𝜓 + 𝑔𝑧) + 𝐷𝑣ℎ𝑐𝑣

′
𝑀𝑊

𝑅𝑇
𝛻𝜓 + 𝐷𝑣ℎ𝑠𝛻𝑇) (15) 

(𝐶𝑚𝜌𝑚𝜒𝑚 + 𝐶𝑤𝜌𝑤𝜃(𝜓))
𝜕𝑇

𝜕𝑡
=  𝛻 ∙ (𝜆𝛻𝑇 + (𝑇𝐶𝑤)𝐾(𝜓)𝛻(𝜓 + 𝑔𝑧) + (𝐿 + 𝑇𝐶𝑤)𝐷𝑣ℎ𝑐𝑣

′
𝑀𝑊

𝑅𝑇
𝛻𝜓  

+(𝐿 + 𝑇𝐶𝑤)𝐷𝑣ℎ𝑠𝛻𝑇) (16) 

2.4 Boundary conditions 

The soil surface conditions and tank temperature determine the moisture and thermal energy flow at the 

interface of soil and atmosphere. The atmospheric temperature 𝑇0 is set as 20
 o

C, and 7
 o

C is constant at the 

lower surface underground. The tank temperature 𝑇𝑡𝑎𝑛𝑘 is fixed as 35
 o

C . Relative humidities of the top and 

bottom surface of underground domain are set as 0.99 and 0.999, respectively. The convective heat transfer 

coefficient between upper boundary and atmosphere ℎ1 is 0.5 W/(m2∙K), for that between tank and soil ℎ2 300 

W/(m2∙K).   

The boundary conditions take the following form regarding other interfaces. 

𝜕𝑇

𝜕𝑥
(𝑊, 𝑦) = 0 (17) 

𝜕𝑇

𝜕𝑦
(𝑥, −𝐻) = 0 (18) 

𝜕𝑇

𝜕𝑦
(𝑥,0) = ℎ1(𝑇 − 𝑇0) (19) 

For the tank: 

868



𝜕𝑇

𝜕𝑟
(𝑟𝑖 , 𝛼) = ℎ2(𝑇 − 𝑇𝑡𝑎𝑛𝑘) (20) 

3. Validation of the non-isolated tank model 

In the coupled simulation model, parameters for calculating the temperature of non-insulated underground 

thermal storage tank and soil around are given in Table 1. As shown in Part 2, the nodes numbering used in the 

calculations are marked in Figure 3. Actual annual weather data from December 2020 to December 2021 was 

inputted into the model to obtain more precise results. The model was operated for one-year totally. Schematic 

of the UHSU experimental measurements is arranged in Figure 4. 

The initial soil temperature around the tank, initial relative humidity of top and bottom layers underground, initial 

top and bottom tank temperature, along with weather data were collected in the first days of selected months to 

verify the coupled model. Two series of validation were conducted: verification of the top tank temperature and 

tank temperature. 

 

Figure 3: Typical nodes in MATLAB simulation for modelling of the UHSU 

 

Figure 4: Schematic of temperature sensors positions in the UHSU experiment 

Table 1: Parameters for UHSU model in MATLAB 

Parameter  Description Parameter Value 

𝑟𝑜 Outer radius of the tank 1.5 m 

𝑊 Ground width 12 m 

𝐻1 Tank burial depth 5 m 

𝐻2 Underground depth 15 m 

ℎ1 Air-to-insulated ground heat transfer coefficient (external surface) 0.5 W/(m2∙K) 

ℎ2 Water-ground heat transfer coefficient in the non-insulated tank 300 W/(m2∙K) 

𝜌𝑤 Water density 1,000 kg/m3 

𝑐𝑤 Water specific heat 4,180 J/(kg∙K) 

𝑐 Ground specific heat 1,000 J/(kg∙K) 

𝜌 Ground density 1,200 kg/m3 

𝑘 Ground thermal conductivity 1.2 W/(m∙K) 

𝐴𝑃𝑉 The surface area of photovoltaic PV panels 42.66 m2 

𝜂𝑃𝑉 Thermal efficiency of PV panels 0.45 

869



For the verification of the top tank temperature, temperature recordings of node 2 in model (Figure 3) and node 

14, 18 and 22 from measurements (Figure 4) were collected to compare. The results are shown in Table 2. 

For the verification of average tank temperature, temperature recordings of node 1, 2, 3, 11 and 12 in model 

(Figure 3) and the average temperature of top and bottom layer of tank in truth measurements were utilised to 

validate. The results are displayed in Table 3. 

Table 2: Validation results of top tank temperature 

Date Temperature in model (°C) Temperature in measurement (°C) Relative error (%) 

2021-05-01 15.532 15.133 2.637 

2021-09-01 28.227 28.300 0.258 

2021-10-01 27.013 26.133 3.367 

2021-11-01 20.886 20.200 3.396 

Table 3: Validation results of average tank temperature 

Date Temperature in model (°C) Temperature in measurement (°C) Relative error (%) 

2021-07-01 34.384 33.550 2.485 

2021-08-01 34.535 33.600 2.782 

2021-09-01 27.592 29.300 5.829 

2021-10-01 26.591 27.550 3.479 

From the results above, one can observe that relative error between model and experiment of the top tank 

temperature varies from 0.258 % (September) to 3.396 % (November) whilst it is in the range of 2.485 % (July) 

to 5.829 % (September) in the case of tank temperature. A possibility account for temperature discrepancy is 

that the simulation process considers little irreversible loss, which results in a fully heat transfer between the 

tank and soil in the underground occurred. Besides, the temperature in the simulation model shows the same 

trend with it in the actual measurement. 

4. Conclusions 

This paper presents the validation process of a coupled model between moisture and heat transfer of RESHeat 

underground heat storage unit. The ordinary differential equations are formulated to obtain the thermal energy 

balance of the underground energy storage system. Actual measurements including non-insulated thermal 

storage tank temperature, soil temperature, relative humidity and weather data have been initialised the 

proposed model. Two types of tank temperatures are analysed between the simulation and experiments. 

Results show that relative error of model and measurement varies from 0.258 % to 3.396 % (for top tank 

temperature) and from 2.485 % to 5.829 % (for average tank temperature). 

References 

Daneshazarian R., Bayomy A.M., Dworkin S.B., 2022, NanoPCM based thermal energy storage system for a 

residential building, Energy Conversion and Management, 254, 115208. 

Kroener E., Vallati A., Bittelli M., 2014, Numerical simulation of coupled heat, liquid water and water vapor in 

soils for heat dissipation of underground electrical power cables, Applied Thermal Engineering, 70, 510-523. 

Li K., Jiang Y., Jia L., Zhao G., Kaita M., 2023, Modeling study on nearly-zero carbon cooling in single houses 

with underground cold storage, Journal of Energy Resources Technology, 145, 012103. 

Ocłoń P., Chin H.H., Kozak-Jagieła E., Taler J., Ścisłowicz F., Czamara M., 2023, Photovoltaic–thermal waste 

heat integration with underground thermal energy storage and heat pump systems, Handbook of Process 

Integration (PI), Woodhead Publishing Limited, Cambridge, UK, 1017-1042. 

Ocłoń P., Cisek P., Kozak-Jagieła E., Taler J., Taler D., 2020, Modeling and experimental validation and thermal 

performance assessment of a sun-tracked and cooled PVT system under low solar irradiation, Energy 

Conversion and Management, 222, 113289. 

Ocłoń P., Ławryńczuk M., Czamara M., 2021, A new solar assisted heat pump system with underground energy 

storage: modelling and optimisation, Energies, 14, 5137. 

Xu Y., Zeng Z., 2022, Experimental and numerical investigation on the effect of heat and moisture coupling 

migration of unsaturated lateritic clay for the soil thermal storage system, Energy and Buildings, 276,112499. 

Zhang S., Ocłoń P., Klemeš J.J., Michorczyk P., Pielichowska K., Pielichowski K., 2022, Renewable energy 

systems for building heating, cooling and electricity production with thermal energy storage, Renewable and 

Sustainable Energy Reviews, 165, 112560. 

870