Numerical Investigation of Harvesting Solar Energy and Anti-Icing Road Surfaces Using a Hydronic Heating Pavement and Borehole Thermal Energy Storage

Abstract

Hydronic Heating Pavement (HHP) is an environmentally friendly method for anti-icing the roads. The HHP system harvests solar energy during summer, stores it in a Seasonal Thermal Energy Storage (STES) and releases the stored energy for anti-icing the road surface during winter. The aims of this study are to investigate: (i) the feasibility of HHP system with low fluid temperature for harvesting solar energy and anti-icing the road surface; and (ii) the long-term operation of the STES. In this study, a Borehole Thermal Energy Storage (BTES) is considered to be the STES. The HHP system and the BTES are decoupled from each other and their performances are investigated separately. A hybrid 3D numerical simulation model is developed to analyze the operation of the HHP system. Moreover, a 3D numerical simulation model is made to calculate the temperature evolution at the borehole walls of the BTES. The climate data are obtained from Östersund, a city in the middle of Sweden with long and cold winter periods. Considering the HHP system with the inlet fluid temperature of 4 °C, the road area of 50 m × 3.5 m as well as the BTES with 20 boreholes and 200 m depth, the result showed that the harvested solar energy during summer is $352.1\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$, the required energy for anti-icing the road surface is $81.2\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ and the average temperature variation at the borehole walls after 50 years is +0.5 °C. Installing the HHP system in the road leads to a 1725 h shorter remaining number of hours of slippery condition on the road surface during winter and a 5.1 °C lower temperature on the road surface during summer, compared to a road without the HHP system.

1. Introduction

Approximately 50% of all annual fatal accidents of passenger cars in Sweden occur in slippery conditions [1]. In rural roads of northern Sweden, the proportion of fatal accidents associated with snow/ice covered roads is about 90% [2]. Hence, the effective winter maintenance is a vital service to ensure the safety and accessibility of roads. A traditional method for de-icing and anti-icing the slippery conditions on the road surface is to spread out salt and sand [3]. De-icing is an action to remove the available ice from the road surface, while anti-icing is an action to prevent ice formation on the road surface [4]. Sanding and salting can cause the pollution of the surface and water ground [5] as well as the corrosion of road infrastructures and vehicles [6]. The consumptions of sand, used for winter road maintenance in Scandinavian countries, is over 600,000 ton and the consumptions of salt is over 1,700,000 ton [7]. Sometimes, sand and salt agents are distributed on the road surface after the presence of slippery condition [8]. Delay in the mitigation of the slippery conditions will result in the occurrence of severe accidents [9]. The number of traffic accident on a snow/ice covered road in Sweden is approximately five times more than that on a dry road surface [10]. Furthermore, absorption of solar radiation during summer months can increase up the temperature of the road surface, consisting of asphalt pavement, up to 70 °C [11]. The high temperature of the road surface usually degrades the durability of the asphalt pavement by accelerating the thermal oxidation in the asphalt pavement [12] and plastic deformations under traffic loads [13].

A renewable alternative method to overcome the above-mentioned problems is using Hydronic Heating Pavement (HHP) [14]. The HHP system includes the embedded pipe in the road, connected to a Seasonal Thermal Energy Storage (STES). The HHP system harvests solar energy during summer and stores it in the storage. During winter, the stored energy, as the only source of energy, is injected into the embedded pipes to warm up the road surface.

The earliest hydronic heating system was constructed in Oregon, USA in 1948. The system worked based on geothermal hot water and was in operation for 50 years before it collapsed due to the external corrosion of metal pipes [15]. The Solar Energy Pilot Project (SERSO), constructed on a bridge in Därligen in Switzerland in 1994 [16], and the hydronic heating system of Göteborgsbacken, constructed on a road section close to the city of Jönköping in Sweden in 2007 [17], are two successful examples for the heating systems in roads. The idea of both projects were to keep the surface temperature above 0 °C to prevent ice formation and freezing the compacted snow on the bridge/road surface [18]. The SERSO project, which consists of 91 boreholes with 65 m deep, annually run for less than 1000 h in summer to harvest solar energy and another 1000 h in winter to mitigate the slippery conditions on the bridge surface [19]. Göteborgsbacken project, which is connected to a district heating system, can provide 125–250 $\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ energy for anti-icing the road surface [17].

To analyze the operation of the HHP system, different numerical simulation models were developed. Pahud [20] made a one-dimensional (1D) numerical model to simulate the heating system, installed in a bridge. This heating system was controlled using the air temperature. The system was turned on when the air temperature was below 4 °C and turned off when the air temperature was below −8 °C. The reason for turning off the system was the low likelihood of the ice formation on the road surface due to the low humidity content of air for the air temperature below −8 °C [19]. Moreover, Abbasi [21] developed a 1D numerical model of the HHP system in order to investigate the anti-icing operation for two modes: (i) keeping the road surface temperature above 0 °C and (ii) keeping the road surface temperature above the dew-point temperature when the surface temperature is below 0 °C. Abbasi found that applying the second mode resulted in ten times less annual required energy for anti-icing the road surfaces.

Borehole Thermal Energy Storage (BTES) is one of the most well-known technologies for application of the HHP systems [15]. BTES uses soil/rock as the heat storage medium and consists of vertically drilled holes with a depths of 30 m to 200 m [22]. Johnsson [23] investigated the feasibility of the coupled HHP system to the BTES for ten different cities in Sweden. By assuming that 30% of the harvested solar energy during summer can be re-used in winter for heating the road surface, Johnsson [23] showed that the mean value of the harvested energy is 131 $\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ with a standard deviation of 7 $\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ and the mean value of the required energy is 119 $\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ with a standard deviation of 73 $\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$. The high standard deviation for the required energy indicated that the local climate conditions can cause rather larger variations in the energy demand [23].

There are certain building energy simulation platforms such as “BRIDGESIM” build on TRNSYS [20] and “Geothermal smart bridge” build on HVACSIM+ [24] to simulate the whole system including HHP system, heat pumps, short term storage and BTES [22]. However, there is a general inadequacy and limitation for these platforms [25]. For example, the HHP system was considered to be either a simplified 1D model or a 2D model and the influence of wind, precipitation and latent heat on the road surface were neglected [23]. In this study, instead, a hybrid three-dimensional (3D) numerical simulation model is developed to simulate the transient operation of the HHP system. The hybrid 3D model simultaneously calculates the transient heat, flowing out from the pipes based on the Finite Element Method (FEM) and the fluid temperature decline along the pipe. Furthermore, a 3D numerical simulation model of the BTES is created to find out the temperature evolution in the borehole walls on long-term operation. It is worth mentioning that, the authors of this paper already simulated the harvesting and anti-icing operations of the HHP system [14,26]. However, the studies were based on a 2D model where the calculation of the fluid temperature decline along the pipe was not considered.

The aims of this study are: (i) to investigate the feasibility of the HHP system with low temperatures for harvesting solar energy during summer and anti-icing the road surface during winter as well as (ii) to investigate the feasibility of the BTES for injection/extraction heat from the ground on long-term operation.

In practice, heat pumps are required for coupling the HHP system to the BTES. The heat sink and source of BTES can be affected by the heating and cooling operation of heat pump system. However, in this study, only the basic concept of using the HHP system and BTES for harvesting solar energy and anti-icing the road surface are investigated and the design and the performance of the heat pumps are not studied in detail. The HHP system and the BTES are decoupled from each other and their performances are investigated separately. The scheme of the HHP system, heat pump and BTES are shown in Figure 1.

The numerical simulation models were made in COMSOL Multiphysics 5.3. The simulation was run in a computer with 16 GB RAM and a processor of Intel (R) Core (TM) i7-4600K CPU @ 2.1 GHz. To determine if the mesh size in the FEM is small enough to get accurate results, the numerical simulation models were run with three different mesh sizes with the maximum element size of 0.103 m, 0.038 m, and 0.012 m. The relative error between the results associated with temperature change on the road surface and the outlet fluid temperature for the three different mesh sizes were less than 1%. In this study, the mesh size with the maximum element size of 0.069 m and the minimum element size of $3.07\times {10}^{-4} \mathrm{m}$ was used to make the numerical simulation models. It took about 60 min to run the hybrid 3D numerical simulation model of the HHP system and about 250 min to run the 3D numerical simulation model of the BTES. Furthermore, in order to check the accuracy of the hybrid 3D numerical simulation model, the numerical simulation model was validated using the results obtained from the literature [27]. The validation results are presented in Section 3.2. For more details of the validation of the hybrid 3D numerical simulation model, the reader is referred to [28,29,30].

2. Mass and Heat Balance

This section describes the mass and heat balances on the road surface. The equations are based on [26,31]. In this study, it is assumed that all snowfalls are immediately plowed away from the road surfaces, the rainfall and the water produced due to the melted ice/snow are drained well and the amount of sublimation and deposition are negligible. The mass balance of the water on the road surface is only based on condensation and evaporation. If $m^{″}{}_{moisture}(\mathrm{kg}/{\mathrm{m}}^2)$ is the mass balance of the water on the road surface, ${\dot{m}}^{″}{}_{con}(\mathrm{kg}/({\mathrm{m}}^2·\mathrm{s}))$ is the condensation rate and ${\dot{m}}^{″}{}_{evp}(\mathrm{kg}/({\mathrm{m}}^2·\mathrm{s}))$ is the evaporation rate, the mass balance of water on the road surface is obtained as:

$$\frac{d{m}_{moisture}^{″}}{dt}={\dot{m}}_{con}^{″}-{\dot{m}}_{evp}^{″}$$

(1)

It is important to note that the mass balance is considered to be a control point for the heat balance; i.e., the required heat for the evaporation has to be stopped when the amount of evaporated water from the road surface is equal to the amount of the existing water on the road surface due to the condensation.

Furthermore, the heat balance on the road surface consists of seven heat fluxes as Equations (2)–(8). It should be noted that since all snowfalls are assumed to be removed from the road surface, then the heat flux associated with the latent heat of the snow melting is not taken into account. However, falling snow, before the snow removal is started, affects the heat balance of the road surface. Hence, the sensible heat flux of snow, $q_{snow}$ (W/m²), related to the heat capacity of snow and temperature change is considered in the heat balance:

• conductive heat from ground and pipes:

  $$q_{cond}=-\lambda ·\nabla T$$

  (2)
• convective heat flow from the ambient air

  $$q_{conv}=h_c·(T_{ambient}-T_{surface})$$

  (3)
• sensible heat from rain

  $$q_{rain}={\dot{m}}_{rain}·c_{p-water}·(T_{ambient}-T_{surface})$$

  (4)
• sensible heat from snow

  $$q_{snow}={\dot{m}}_{snow}·c_{p-snow}·(T_{ambient}-T_{surface})$$

  (5)
• long-wave radiation

  $$q_{lw}=\epsilon ·\sigma ·(T_{sky}^4-T_{surface}^4)$$

  (6)
• short-wave radiation

  $$q_{sw}=\alpha ·I$$

  (7)
• latent heat of evaporation and condensation

  $$q_{evp/con}=h_e·\beta ·(v_{ambient}-v_s)$$

  (8)

  where $\lambda$ (W/(m·K)) is the thermal conductivity of the road materials, $T$ (K) is the temperature, $T_{ambient}(\mathrm{K})$ is the ambient air temperature, $T_{surface}(\mathrm{K})$ is the road surface temperature, $h_c(\mathrm{W}/({\mathrm{m}}^2·\mathrm{K}))$ is the convective heat transfer coefficient, $\epsilon$ (-) is the emissivity of the surface, $\sigma (\mathrm{W}/({\mathrm{m}}^2·{\mathrm{K}}^4))$ is the Stefan-Boltzmann constant, $T_{sky}(\mathrm{K})$ is the sky temperature, $\alpha$ (-) is the solar absorptivity of the surface, $I(\mathrm{W}/{\mathrm{m}}^2)$ is the solar irradiation, $h_e(\mathrm{J}/\mathrm{kg})$ is the latent heat of evaporation of water, $\beta (\mathrm{m}/\mathrm{s})$ is the moisture transfer coefficient, $v_{ambient}(\mathrm{kg}/{\mathrm{m}}^3)$ is the humidity by the volume of the ambient air, $v_s(\mathrm{kg}/{\mathrm{m}}^3)$ is the humidity by the volume of the saturated air at the surface temperature and $\dot{m}(\mathrm{kg}/({\mathrm{m}}^2·\mathrm{s}))$ is the mass rate of snow/rain per square meter of the surface. For more details about the calculation of the different parameters, the reader is referred to [14].

The heat balance of the road surface is calculated as:

$$q_{cond}+q_{conv}+q_{rain}+q_{snow}+q_{lw}+q_{sw}+q_{evp/con}=0$$

(9)

An illustration of the different heat fluxes on the road surface is shown in Figure 2.

In Figure 2, the directions of heat fluxes are toward the road surface, regardless of their sign. In the rest of the work, the positive sign (+) is used for the heat flux which is given to the surface and the negative sign (−) is used for the heat flux which is taken out from the surface.

3. Numerical Simulation Model

This section presents: (i) development of the hybrid 3D numerical simulation model, (ii) validation of the hybrid 3D model, (iii) boundary conditions of the numerical simulation model, (iv) criteria for anti-icing the road surface and (v) criteria for harvesting solar energy.

3.1. Hybrid 3D Numerical Simulation Model

Hybrid 3D numerical simulation model of the HHP system is developed based on the study done by Karlsson [32], which was related to the thermal modeling of water-based floor heating system for buildings. This section consists of three sub-sections, namely: (i) the spatial discretization of the HHP system, (ii) the heat transfer process between the fluid and the embedding materials and (iii) the coupling of adjacent sections along the pipe. The equations in this section are extracted from [31,32].

3.1.1. Spatial Discretization of the HHP System

The sketch of a hybrid 3D model of the HHP system is shown in Figure 3a). As can be seen, the road surface is subdivided into a finite number of smaller sub-surfaces. In addition, a spatial grid is produced along the embedded pipe. The size of the gird is determined by the distance between pipes. The length of a pipe is also subdivided into the smaller sections, resulting in a sequence of pipe sections. For each subsurface/pipe section, the 3D structure of the road is represented by two 2D vertical sections which are perpendicular to pipe direction, see Figure 3b. Each 2D vertical section includes a road structure and an embedded pipe. All 2D sections are serially connected through the convective heat transfer in the fluid, circulating along the pipes. The top boundary of each 2D section is bound to the ambient heat transfers and the bottom boundary is bound to a constant temperature (details are explained in Section 3). The left and right boundaries are considered to be adiabatic.

3.1.2. Heat Transfer Process between the Fluid and Embedding Materials

This section illuminates the essential heat transfer process between the fluid and the embedding materials. A 2D view of the road section including the embedded pipe is presented in Figure 4a. Quadrilateral elements are used to create the mesh on the road domain. As can be seen in Figure 4b, the positions of nodes are indicated by index i in the horizontal direction and by index j in the vertical direction. The total thermal resistance between the fluid and the embedding materials is represented by $R_{eq-pipe} (\mathrm{m}·\mathrm{K}/\mathrm{W})$. $R_{eq-pipe}$ consists of three serially coupled thermal resistance, namely: (i) the surface heat transfer resistance at inner pipe surface, $R_{PWS} (\mathrm{m}·\mathrm{K}/\mathrm{W})$, (ii) the pipe resistance, $R_p (\mathrm{m}·\mathrm{K}/\mathrm{W})$, and (iii) the thermal resistance between the outer pipe wall and an additional annulus, $R_{i.j} (\mathrm{m}·\mathrm{K}/\mathrm{W})$. The value of $R_{PWS}$ is calculated using Reynolds, Prandtl and Nusselt dimensionless numbers, see Equations (10) to (13).

$$\mathrm{Re}=\frac{2·{\rho }_f·v_f·r_{inner}}{{\mu }_f}$$

(10)

$$\mathrm{Pr}=\frac{{\mu }_f·c_{p,f}}{{\lambda }_f}$$

(11)

$$\mathrm{Nu}=\{\begin{array}{c}4 \mathrm{if} \mathrm{Re}\le 2300\\ 0.023·{\mathrm{Re}}^{0.8}·{\mathrm{Pr}}^{1/3} \mathrm{if} \mathrm{Re}>2300\end{array}$$

(12)

$$R_{PWS}=\frac{1}{\pi ·{\lambda }_f·\mathrm{Nu}}$$

(13)

where Re (-) is the Reynolds number, Pr (-) is the Prandtl number, Nu (-) is the Nusselt number, ${\rho }_f (\mathrm{kg}/{\mathrm{m}}^3)$ is the density of fluid, $v_f (\mathrm{m}/\mathrm{s})$ is the average velocity of fluid through a cross section of pipe, $r_{inner} (\mathrm{m})$ is the inner radius of the pipe, $u_f (\mathrm{kg}/(\mathrm{m}·\mathrm{s}))$ is the dynamic viscosity of fluid, $c_{p,f} (\mathrm{J}/(\mathrm{kg}·\mathrm{K}))$ is the specific heat capacity of fluid and ${\lambda }_f (\mathrm{W}/(\mathrm{m}·\mathrm{K}))$ is the thermal conductivity of fluid. Furthermore, $R_p$ is calculated as:

$$R_p=\frac{\ln (\frac{r_{outer}}{r_{inner}})}{2·\pi ·{\lambda }_{pipe}}$$

(14)

where $r_{inner} (\mathrm{m})$ and $r_{outer} (\mathrm{m})$ are respectively the inner and outer radii of the pipe and ${\lambda }_{pipe} (\mathrm{W}/(\mathrm{m}·\mathrm{K}))$ is the thermal conductivity of the pipe material. Finally, $R_{i,j}$ is calculated as:

$$R_{i,j}=\frac{\ln (\frac{r_{i,j}}{r_{outer}})}{2·\pi ·{\lambda }_{i,j}}$$

(15)

where $r_{i,j} (\mathrm{m})$ is the radius of the additional annulus surrounding the pipe and ${\lambda }_{i,j} (\mathrm{W}/(\mathrm{m}·\mathrm{K}))$ is the thermal conductivity of the surrounding materials. It should be noted that $r_{i,j}>r_{outer}$. The total thermal resistance between the fluid and the embedding materials, $R_{eq-pipe}$, is obtained from Equation (16). Moreover, the equivalent temperature surrounding the pipe, $T_{eq-pipe} (\mathrm{K})$, corresponds to the average ambient temperature of the embedding materials close to the pipe segment is calculated as Equation (17).

$$R_{eq-pipe}=R_p+R_{i,j}+R_{PWS}$$

(16)

$$T_{eq-pipe}=T_f-q_{i,j}·R_{eq-pipe}$$

(17)

where $T_f (\mathrm{K})$ is the temperature of fluid, circulating through the pipe and $q_{i,j} (\mathrm{W}/{\mathrm{m}}^2)$ is the heat flow between the fluid and the embedding materials.

3.1.3. Coupling Adjacent Sections along the Pipe

The longitudinal fluid temperature distribution and the transversal heat transfer process along the pipe are calculated using a steady state solution as Equation (18).

$$v_f·\pi ·r_{inner}{}^2·{\rho }_f·c_{p,f}·\frac{\partial T_f(x)}{\partial x}+\frac{T_f(x)-T_{eq-pipe}(x)}{R_{eq-pipe}}=0$$

(18)

where x (m) is the longitudinal length coordinate along the fluid motion, $T_f(x)$ and $T_{eq-pipe}(x)$ are respectively the fluid temperature and the equivalent temperature surrounding the pipe at the given point of x.

The distribution of the longitudinal fluid temperature along the pipe is calculated considering the following assumptions:

• The transversal heat transfer flow for each 2D section of the HHP system is calculated independently of the adjacent section.
• The fluid temperature distribution is calculated using a local quasi steady state condition. The condition is valid for each time step.
• The fluid temperature distribution is calculated using a step-wise sequence solution for all 2D sections of the HHP systems.
• The conductive heat transfer of the fluid inside the pipe is neglected since it is small in comparison with the convective heat transfer along the pipe.

As shown in Figure 5, the total length of the pipe is subdivided into a finite number of sections. The sections are numbered with the index of k, starts from 1 and continues along the fluid motion. The length of pipe sections is labeled as $L_1$, $L_2$, … $L_n$.

For the pipe section between the labels of n and n + 1, the section length is $L_n$, the inlet temperature of the fluid at time t is $T_{f,n}^t$ and the outlet temperature is $T_{f,n+1}^t$. Considering the steady state longitudinal fluid temperature distribution along the pipe, the analytical solution to obtain the outlet fluid temperature is as:

$$T_{f,n+1}^t=T_{eq-pipe,n}^t+(T_{f,n}^t-T_{eq-pipe,n}^t)·e^{-(L_n/l_n)}$$

(19)

where $l_n (\mathrm{m})$ is the characteristic length related to the interaction between the convective heat transfer along the pipe and the transversal heat flow. The value of $l_n$ is calculated as:

$$l_n=R_{eq-pipe}·v_f·\pi ·r_{inner}{}^2·{\rho }_f·c_f$$

(20)

In this study, the inlet temperature of fluid, $T_{f,1}^t$ (K), is a known parameter. By applying Equations (19) and (20), it is possible to calculate the $T_{f,2}^t$ (K). The value of $T_{f,2}^t$ will be used as the inlet temperature for the second pipe section. By considering the outlet temperature from one section as the inlet temperature for the next section, the outlet temperature for whole length of the pipe will be equal to the outlet temperature at the last pipe section.

3.2. Validation of the Numerical Model of the HHP System

The hybrid 3D numerical simulation model, developed in this study, is validated by the results of another numerical simulation model from the literature [27] which was about investigation of heat-collecting properties of asphalt pavement as a solar collector by a three-dimensional unsteady model. The solar collector consisted of a one-layer asphalt pavement and the embedded pipes. The dimension of the solar collector was 4 m × 4 m × 0.7 cm (length × width × thickness). The pipes were in a serpentine configuration. The total length of pipe was 63.15 m. The distance between pipes was 0.25 m (center to center) and the embedded depth of pipes was 60 mm (from surface to center). The fluid was water. The velocity of fluid, circulating inside the pipes, was 0.3 m/s. The inlet temperature of the fluid was 18 °C. The surface boundary was exposed to the solar radiation, convection, and long-wave radiation. The solar radiation had the constant value of 1080 $\mathrm{W}/{\mathrm{m}}^2$. The bottom and the side boundaries of the solar collector were considered to be adiabatic. The schematic view of the solar collector is shown in Figure 6. For more details, the reader is refered to [27].

The thermal properties of materials including thermal conductivity, density and specific heat capacity are given in

Table 1. Thermal properties of the concrete pavement, pipe and fluid (water) [27].

Material	$\mathrm{Thermal}\mathrm{Conductivity}(\mathrm{W}/(\mathrm{m}·\mathrm{K}))$	Density $(\mathrm{k}\mathrm{g}/{\mathrm{m}}^3)$	Specific Heat Capacity $(\mathrm{J}/(\mathrm{k}\mathrm{g}·\mathrm{K}))$
concrete pavement	1.62	2339	896
Pipe	0.42	932	2300
Fluid (water)	0.61	1000	4181

.

The outer and inner diameters of pipe were assumed to be 20 mm and 15 mm. Furthermore, the value of absorptivity and emissivity were set to be 0.85 (-). The valiation were done for two cases:

• Case 1: the ambient air temperature was 25 °C and the initial temperature was 20 °C.
• Case 2: the ambient air temperature was 34 °C and the initial temperature was 30 °C.

The validation was done based on: (i) the temperature of aspahlt pavement at the depth of 2.5 cm and (ii) the outlet temperature of fluid. From the literature [27], the temperature variation was given for 120 min. Hence, to validate the hybrid 3D numerical simulation model, the harvesting process was run for the same time with the time step of 1 min. The validation results are shown in Figure 7.

As can be seen from Figure 7a, the temperature of the asphalt pavement at the depth of 2.5 cm obtained from the literature [27] and the hybrid 3D numerical simulation model in this study are matching well. The absolute temperature difference between the two models for the case 1 is 0.3 °C with the standard deviation of 0.2 °C and that for the case 2 is 0.2 °C with the standard deviation of 0.2 °C. Furthermore, from Figure 7b, the outlet temperature of the fluid, obtained from the literature [27] and the hybrid 3D numerical simulation model in this study have good match with together, so the absolute temperature difference between the two models for the case 1 is 0.3 °C with the standard deviation of 0.2 °C and that for the case 2 is 0.4 °C with the standard deviation of 0.3 °C. Considering the results of validation, the hybrid 3D numerical simulation model is used in the rest of this study to investigate the operation of the HHP system for harvesting solar energy and anti-icing road surface.

3.3. Boundary Conditions of the Numerical Simulation Model

The hybrid 3D numerical simulation model simultaneously calculates the transient heat, flowing out from the pipes based on the FEM and the fluid temperature decline along the pipe. If it is assumed that the equal amount of heat fluxes is incoming and outgoing from the boundaries of A-B and C-D in Figure 2, the region of A-B-C-D can be considered as a symmetrical analysis region. To reduce the computation time, this symmetrical region is selected to create the numerical simulation model of the HHP system.

The boundary of B-C is exposed to the mass and heat balances based on Equations (1)–(8). Two boundaries of A-B and C-D are set to be adiabatic. Furthermore, the boundary of A-D, which is located at the depth of 5 m from the road surface, is set to have a constant temperature of 2.53 °C. This temperature is equal to the annual mean temperature of ambient air in Östersund.

Furthermore, as can be seen from Figure 3b, the hybrid 3D model of the HHP system was represented by four 2D vertical cross sections. For details about selecting the number of cross sections, the reader is referred to [28].

In this study, the climate data are obtained from Östersund, a city in middle of Sweden. The city has long and cold winter period which is interesting to simulate the anti-icing operation of the HHP system. The location of Östersund is shown in Figure 8.

The climate data were generated at a 1 h interval [33] and included: dry-bulb/air temperature (°C), relative humidity (%), wind speed (m/s), dew-point temperature (°C), incoming long-wave radiation (W/m²), short-wave radiation (W/m²) and precipitation (mm/h). The numerical simulation was modeled using the implicit time-stepping method. The time step was 1 h in line with the climate data resolution.

The road structure consists of six different layers. The first three layers are asphalt pavement, the next two layers are crushed aggregates and the last layer is ground. The material type, thickness and the thermal properties of each layer are presented in

Table 2. Materials properties of the road layers.

Material	Thickness (mm)	Thermal Conductivity (W/(m·K))	Density (kg/m3)	Specific Heat Capacity (J/(kg·K))
Wearing layer (ABS16)	40	2.24	2415	848
Binder layer (ABB22)	60	1.44	2577	822
Base layer (AG22)	100	1.51	2582	894
Subbase layer	80	0.7	1700	900
Subgrade layer	1000	0.8	1400	900
Ground	3720	0.6	1300	600

. The thermal properties related to the wearing layer, binder layer and base layer were experimentally measured by authors in a room with a fixed temperature of $22°\mathrm{C}\pm 0.2°\mathrm{C}$ [14,34]. The thermal properties of the subbase layer, subgrade layer and the ground are obtained from literatures [31,35,36,37,38]. The thermal properties of the materials were reported in the temperature range between $20°\mathrm{C}\mathrm{and}25°\mathrm{C}.$

Moreover, the data related to the absorptivity and emissivity of the road surface, the geometrical and thermal properties of the pipe materials as well as the thermal properties of the fluid are presented in

Table 3. Information about the initial simulated HHP system in this study (Fluid is 25% ethylene glycol-water mixture [39]).

Parameter	Value	Unit
Thermal conductivity of pipe materials	0.4	W/(m·K)
Density of pipe materials	925	kg/m3
Specific heat capacity of pipe materials	2300	J/(kg·K)
Outer diameter of the embedded pipes	25	mm
Inner diameter of the embedded pipes	20.4	mm
Diameter of additional annulus ($2·r_{i,j}$)	30.4	mm
Distance between the pipes (center to center)	100	mm
Embedded depth (from center of the pipe to the surface)	87.5	mm
Emissivity of the road surface	0.89	-
Absorptivity of the road surface	0.78	-
Pipe length	50	m
Fluid flow rate	8	l/min
Thermal conductivity of fluid	0.488	W/(m·K)
Specific heat capacity of fluid	3852	J/(kg·K)
Density of fluid	1025	kg/m3
Dynamic viscosity	2.9	mPa·s

. The data associate with the pipe material, made from polyethylene (PEX) and measured in the temperature of 20 °C, are obtained from [37]. The data associate with the thermal properties of fluid, 25% ethylene glycol-water mixture, are extracted from [38] for the temperature range between 25 °C and 30 °C.

3.4. Criteria for Anti-Icing the Road Surface

When the surface temperature, $T_{surface}$ (K), is below the dew-point temperature, $T_{dew}$ (K), the water will be formed on the road surface due to the condensation. In addition, if $T_{surface}$ is below the freezing temperature of water, $T_{freezing}$ (K), the water on the road surface will turn to ice. An idea to avoid the occurrence of slippery conditions is to keep the surface temperature above $T_{dew}$, when $T_{surface}<T_{freezing}$. The mentioned criteria for running the heating system could be written as:

$$\{\begin{array}{c}{T}_{surface}<T_{dew}\\ T_{surface}<T_{freezing}\end{array}$$

(21)

Due to the non-uniform distribution of the heat fluxes on the road surface after turning on the heating system, the slippery conditions are not mitigated simultaneously on all points of the road surface. The points on the road surface which are directly above the heating pipe (points B and F in Figure 2) will have the maximum temperature on the road surface and the shortest number of hours of the slippery conditions, compared to the other points between B and F. On the contrary, the minimum surface temperature and the longest number of hours of the slippery conditions will occur on the surface in the middle between two pipes (point C in Figure 2). Furthermore, the fluid temperature will decline as circulating inside the pipe. Hence, the coldest point on the road surface will be at the outlet section and in the middle between two pipes, see “Control point” in Figure 3b.

Whenever the heating system is turned on, the temperature of the fluid is set to be a constant temperature. Furthermore, whenever the heating system is turned off, the boundary condition at the inner surface of the pipe wall is set to be adiabatic. The annual required energy for anti-icing the road surface using the temperature difference between the inlet and outlet fluids, $E_r$ (kWh/(m²·year)), is calculated as:

$$E_r=\frac{1}{c·L}·{{\displaystyle \int }}_{t=0}^{1year}(T_{f,out}^t-T_{f,in}^t)·v_f·\pi ·r_{inner}{}^2·{\rho }_f·c_{p,f}·dt$$

(22)

where c (m) is the distance between pipes, L (m) is the pipe length, $T_{f,in}^t$ (K) is the inlet temperature of the fluid and $T_{f,out}^t$ (K) is the outlet temperature of the fluid, $v_f (\mathrm{m}/\mathrm{s})$ is the velocity of the fluid, $r_{inner} (\mathrm{m})$ is the inner radius of the pipe,${\rho }_f (\mathrm{kg}/{\mathrm{m}}^3)$ is the density of the fluid and $c_{p,f} (\mathrm{J}/(\mathrm{kg}·\mathrm{K}))$ is the specific heat capacity of the fluid.

Furthermore, the number of hours of the slippery condition on the road surface, $t_{slippery}$ (h), which is the number of hours when the temperature of the road surface is lower than both freezing and dew-point temperatures, is calculated as:

$$t_{slippery}={{\displaystyle \int }}_0^{1year}f·dt(f=1\mathrm{if}\left(\begin{array}{c}{T}_{surface}<T_{dew}\\ T_{surface}<T_{freezing}\end{array}\right)\mathrm{else}f=0)$$

(23)

3.5. Criteria for Harvesting Solar Energy

In this study, the air temperature is considered to be the key criterion to start harvesting solar energy. Considering the climate data in Östersund [33], there will be no risk for ice formation on the road surface if the air temperature is above $4°\mathrm{C}$ [26]. Considering a lower air temperature such as $4°\mathrm{C}$ to start harvesting solar energy can result in heating the road surface rather than harvesting. On the contrary, considering a higher air temperature such as $20°\mathrm{C}$, which occurs only 2% during a year, can result in missing many sunny hours for harvesting solar energy. Based on the climate data in Östersund [33], the months from April to September are considered as the harvesting period [30]. For about 70% of this period, the air temperature is above $6°\mathrm{C}$ and for about 10% of this period the air temperature is above $16°\mathrm{C}$ [33]. In this study, five different air temperatures of $6°\mathrm{C}$, $8°\mathrm{C}$, $10°\mathrm{C}$, $12°\mathrm{C}$ and $16°\mathrm{C}$ are arbitrary selected as a criterion to start harvesting solar energy. If $T_{limit} (°\mathrm{C})$ is one of the mentioned five temperatures, then the criterion to start harvesting solar energy can be written as:

$$T_{air}>T_{limit}$$

(24)

Harvesting solar energy will cause a decrease in the temperature of the road surface. The temperature decreases on the road surface during harvesting period, $T_{decrease}$ $(°\mathrm{C})$, can be calculated as:

$$T_{decrease}=T_{unheated}-T_{harv}$$

(25)

where $T_{unheated} (°\mathrm{C})$ is the surface temperature of an unheated road and $T_{harv} (°\mathrm{C})$ is the surface temperature of the HHP system during harvesting period at “Control point”, see Figure 3b.

Similar to Equation (22), the annual harvested solar energy, $E_h (\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year}))$, can be calculated as:

$$E_h=\frac{{{\displaystyle \int }}_{t=0}^{1year}(T_{out-r}-T_{in-r})·v_f·\pi ·r_{inner}{}^2·{\rho }_f·c_{p,f}·dt}{c_r·L_r} (\mathrm{if} T_{air}>T_{limit})$$

(26)

Whenever the harvesting operation is turned on, the temperature of the fluid is set to be a constant temperature. Furthermore, whenever the harvesting operation is turned off, the boundary condition at the inner surface of the pipe wall is set to be adiabatic.

4. Results Related to Harvesting and Anti-Icing Operation of the HHP System

In the previous experimental studies [40,41], the fluid temperature varied at a range between 25 °C and 50 °C. This range of temperature was used to snow melting, however, is high for anti-icing the road surface. The aim of anti-icing is to keep the $T_{surface}$ higher than (but close to) $T_{dew}$ when $T_{surface}<T_{freezing}$, see Section 3.4. Considering that the maximum air temperature in Östersund is $25°\mathrm{C}$ [33], in this study, seven fluid temperatures of 4 °C, 6 °C, 8 °C, 10 °C, 12 °C, 16 °C and 20 °C are taken into account to investigate the feasibility of the HHP system for harvesting solar energy and anti-icing the road surface. The results related to the annual harvested solar energy and the annual required energy for anti-icing the road surface, the average outlet temperatures of the fluid in the HHP system during harvesting and anti-icing periods, the average temperature reduction on the road surface during the harvesting period and the remaining number of hours of the slippery conditions on the road surface during anti-icing period are shown in Figure 9.

As shown in Figure 9a, an increase in the inlet temperature of the fluid, $T_{in}(°\mathrm{C})$, from 4 °C to 20 °C leads to an increase in the annual required energy for anti-icing the road surface, $E_r$, from $81.2 \mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ to $96.3 \mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$. Higher $T_{in}$ causes higher surface temperature and consequently results in more energy consumption. However, an increase in the $T_{in}$ from 4 °C to 16 °C results in a decrease in the annual harvested solar energy, $E_h$, on average from $271.7 \mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ to $34.7 \mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$. Interesting point is that for $T_{in}=20 °\mathrm{C}$, the HHP system is mainly working as a heater rather than a solar collector, so as $44.4 \mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ of energy, on average, is conducted from the pipes to the road surface during harvesting period.

The harvested solar energy highly depends on the value of $T_{limit}$, see Equation (24). By increasing $T_{limit}$ from $6°\mathrm{C}$ to $16°\mathrm{C}$, the time that the HHP system is turned on for harvesting solar energy decreases from 3281 h to 492 h. It should be noted that considering a higher value of $T_{in}$ and a lower value of $T_{limit}$ can result in heating the road surface during harvesting period, rather than harvesting solar energy. For example, considering $T_{in}=20°\mathrm{C}$ and $T_{limit}=6°\mathrm{C}$ results in consuming $140\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ of energy for heating the road surface during harvesting period. This value is even higher than the required energy for anti-icing the road surface, $E_r=96.3 \mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$, when $T_{in}=20°\mathrm{C}$. To harvest maximum possible solar energy from the road surface, it is essential to use a suitable combination of $T_{in}$ and $T_{limit}$. The combinations of $4 °\mathrm{C}\le T_{in}\le 6 °\mathrm{C}$ and $T_{limit}=6 °\mathrm{C}$, $8 °\mathrm{C}\le T_{in}\le 10 °\mathrm{C}$ and $T_{limit}=8°\mathrm{C}$, $T_{in}=12°\mathrm{C}$ and $T_{limit}=10°\mathrm{C}$, $T_{in}=16°\mathrm{C}$ and $T_{limit}=12°\mathrm{C}$ as well as $T_{in}=20°\mathrm{C}$ and $T_{limit}=16°\mathrm{C}$ lead to the maximum values of $E_h$ associated with each inlet fluid temperature, see Figure 9a. The maximum harvested solar energy during summer is $352.1\mathrm{kWh}/({\mathrm{m}}^2·\mathrm{year})$ which is related to $T_{in}=4 °\mathrm{C}$ and $T_{limit}=6 °\mathrm{C}$. It is worth mentioning that for $4°\mathrm{C}\le T_{in}\le 12°\mathrm{C}$, the harvested solar energy can be higher than the required energy for anti-icing the road surface. However, for $T_{in}\ge 16°\mathrm{C}$, the harvested solar energy is always lower than the required energy for anti-icing the road surface.

Figure 9b illustrates the variation of average outlet temperature of the fluid, $T_{out}$, versus $T_{in}$. During both harvesting and anti-icing periods, an increase in $T_{in}$ will results in an increase in $T_{out}$. However, the absolute difference between the inlet and outlet temperatures,$\left|T_{in}-T_{out}\right|$, during harvesting and heating periods is different. During harvesting period and for a constant value of $T_{limit}$, an increase in $T_{in}$ results in harvesting less solar energy, which in turn leads to a lower value of $\left|T_{in}-T_{out}\right|$. On the contrary, during anti-icing period, an increase in $T_{in}$ results in consuming more energy for heating the road surface, which in turn leads to a higher value of $\left|T_{in}-T_{out}\right|$.

Figure 9c illustrates the average temperature reduction on the road surface, $T_{decrease}$, at “Control point”, see Figure 3b. By increasing $T_{in}$ from 4 °C to 20 °C, the value of $T_{decrease}$ follows a decreasing trend from 6 °C to −0.3 °C, on an average value for all $T_{limit}$. $T_{decrease}<0 °\mathrm{C}$, which occurs for $T_{in}=20°\mathrm{C}$ and $6°\mathrm{C}\le T_{limit}\le 10°\mathrm{C}$, means that the temperature of the road surface is increasing during harvesting period, rather than decreasing. The maximum value of $T_{decrease}$ is 6.82 °C which is related to $T_{in}=4°\mathrm{C}$ and $T_{limit}=16°\mathrm{C}$ as well as the minimum value of $T_{decrease}$ is $-1.51 °\mathrm{C}$ which is related to $T_{in}=20°\mathrm{C}$ and $T_{limit}=6°\mathrm{C}$.

Although an increase in $T_{in}$ causes a decrease in $T_{decrease}$, circulating a higher fluid temperature inside the pipe will lead to a decrease in the number of hours of the slippery conditions on the road surface, $t_{slippery}$. From Figure 9d, by increasing $T_{in}$ from 4 °C to 20 °C, the value of $t_{slippery}$ reduces from 284 h to 93 h. It should be noted that the decrease of $t_{slippery}$ follows a convex curve; i.e., as the value of $T_{in}$ increases the number of slippery hours approximates to the preceding value of $t_{slippery}$. It is worth mentioning that using Equation (1) and considering the density of ice as $917 \mathrm{kg}/{\mathrm{m}}^3$ [31], the maximum height of remaining ice on the road surface decreases from 1.48 mm to 0.98 mm as the value of $T_{in}$ increases from 4 °C to 20 °C. The maximum height of the ice on the road surface occurs on the 21st of December at 00:00.

5. Long-Term Operation of the BTES

In this study, it is assumed that harvested solar energy during summer is injected to the BTES and the required energy for anti-icing the road surface during winter is extracted from the BTES. From Figure 9a, it can be seen that for $T_{in}\le 12°\mathrm{C}$, the annual harvested solar energy is higher than the annual required energy for anti-icing the road surface. However, considering the thermal interference between the heat-carrier fluid in the BTES and the undisturbed surrounding ground, the injected heat to the BTES can freely transfer to the surrounding ground [42]. The heat losses from the BTES to the surrounding ground can influence the operation of the BTES in the long term [43]. In this study, to understand the long-term operation of the BTES, the temperature evolutions at the borehole walls are examined over a 50-year period using a 3D numerical simulation model.

To obtain the temperature evolution, the conducted heat flux from the pipes to the road surface during the harvesting and anti-icing periods, $q_p(\mathrm{W}/{\mathrm{m}}^2)$, is extracted from the hybrid 3D numerical simulation model in Section 4. The heat flux of $q_p$ has hourly values in line with the climate data and includes: the harvesting and anti-icing periods as well as the period during which the harvesting/anti-icing is off, $q_p=0 \mathrm{W}/{\mathrm{m}}^2$. The hourly heat fluxes during harvesting and anti-icing periods for different inlet temperatures of the fluids, corresponding to the maximum harvested solar energy in Section 4, are shown in Figure 10.

In Figure 10, the positive values of $q_p$ are related to the harvesting period and the negative values are related to the anti-icing period. The harvesting period is from April to September and the anti-icing period is from October to March. Furthermore, ${\overline{q}}_{p\_year} (\mathrm{W}/{\mathrm{m}}^2)$ is the annual average of heat fluxes at the road surface. As can be seen, the value of ${\overline{q}}_{p\_year}$ is positive and varying from $30.9 \mathrm{W}/{\mathrm{m}}^2$ for $T_{in}=4 °\mathrm{C}$ and $T_{limit}=6 °\mathrm{C}$ to $5.2 \mathrm{W}/{\mathrm{m}}^2$ for $T_{in}=12 °\mathrm{C}$ and $T_{limit}=10 °\mathrm{C}$.

In this study, it is assumed that the harvesting process is done using passive heat exchange, due to the harvested solar heat can transfer from the high temperature road to the low temperature ground. Hence, the amount of the injected heat to the BTES during harvesting process is considered to be equal to the amount of the harvested heat from the road surface. However, it is assumed that the anti-icing process is done using a heat pump to ensure the required inlet temperature to the HHP system and preventing ice formation on the road surface. The extracted heat from the BTES during anti-icing process can be calculated by the difference between the required heat by the HHP system and the produced heat by the heat pump. If ${\mathrm{COP}}_{\mathrm{h}}$ is the coefficient of performance of the heat pump during anti-icing period, the heat flux adjustment factor for anti-icing, $H_f$, can be defined as [44]:

$$H_f=\frac{{\mathrm{COP}}_{\mathrm{h}}-1}{{\mathrm{COP}}_{\mathrm{h}}}$$

(27)

By multiplying the required heat flux of the HHP system by $H_f$, the extracted heat from the BTES can be estimated during the anti-icing process. In this study, for simplicity of the simulation, three constant ${\mathrm{COP}}_{\mathrm{h}}$s of 3, 4 and 5 are used [45]. It is worth mentioning that using a constant ${\mathrm{COP}}_{\mathrm{h}}$ may affect some errors in the hourly simulation results, so by a $1 °\mathrm{C}$ decrease in the entering fluid temperature to the heat pump, the value of ${\mathrm{COP}}_{\mathrm{h}}$ will decrease approximately 0.05 [46]. However, the annual average temperature change in the ground and the thermal interaction of borehole will be unaffected [44].

The injected/extracted heat rate for each borehole, $q_s(\mathrm{W}/\mathrm{m})$, can be obtained using the harvested/anti-icing heat fluxes as:

$$q_s=\frac{W_r·L_r}{d_{total}}·q_p\times f \left(\{\begin{array}{c}f=1 \mathrm{during} \mathrm{harvesting}\\ f=H_f \mathrm{during} \mathrm{anti}-\mathrm{icing}\end{array}\right)$$

(28)

where $L_r(\mathrm{m})$ is the length of the HHP system, $W_r(\mathrm{m})$ is the width of the HHP system and $d_{total}(\mathrm{m})$ is the total required depth of borehole. In this study, $L_r=50\mathrm{m}$ and $W_r=3.5\mathrm{m}$. The value of $d_{total}$ can be determined using the following equation as [47]:

$$d_{total}=\frac{q_h·R_b+q_y·R_{10y}+q_m·R_{1m}+q_h·R_{6h}}{T_m-(T_g+T_p)}(d_{total}=N\times d_{bore})$$

(29)

where N (-) is the number of boreholes, $d_{bore}(\mathrm{m})$ is the depth of each borehole, $q_h(\mathrm{W})$ is the maximum hourly ground peak load, $q_m(\mathrm{W})$ is the average monthly ground load during the month when $q_h$ occurs, $q_y(\mathrm{W})$ is the net heat load of the ground and $R_b(\mathrm{m}·\mathrm{K}/\mathrm{W})$ is the effective borehole thermal resistance. $R_{10y}(\mathrm{m}·\mathrm{K}/\mathrm{W})$, $R_{1m}(\mathrm{m}·\mathrm{K}/\mathrm{W})$ and $R_{6h}(\mathrm{m}·\mathrm{K}/\mathrm{W})$ represent, respectively, the effective ground thermal resistance of yearly, monthly, and hourly ground load components. $T_m(°\mathrm{C})$ is the mean fluid temperature in the borehole, $T_g(°\mathrm{C})$ is the undisturbed ground temperature and $T_p(°\mathrm{C})$ is the temperature penalty. For more details of the parameters, the reader is referred to [47,48].

Considering $R_b=0.1 \mathrm{m}·\mathrm{K}/\mathrm{W}$ and using the average heat loads related to the different fluid temperatures from $4 °\mathrm{C}$ to $12 °\mathrm{C}$ and three different ${\mathrm{COP}}_{\mathrm{h}}$s of 3, 4 and 5, the value of $d_{total}$ is calculated as 3898 m. In this simulation, the depth of each borehole, $d_{bore}$, is considered to be 200 m and it is assumed that the BTES consists of 20 boreholes. Furthermore, it is assumed that the boreholes are located far from each other and there is no thermal interaction among them, so the total heat loads of the BTES are equally divided among all boreholes.

The scheme of the BTES with a single borehole are shown in Figure 11.

In Figure 11, the domain of A-B-C-D-Á-B’-Ć-D’ presents a ground volume which has a rectangular cross section area of 200 m × 200 m. The depth of ground volume is truncated to five time the borehole depth, $5\times d_{bore}$. In Figure 11, the four domains of A-E-O-H-Á-É-Ó-H’, B-E-O-F-B’-É-Ó-F’, C-F-O-G-Ć-F’-Ó-Ǵ, D-G-O-H-D’-Ǵ-Ó-H’ are symmetrical. Hence, to reduce the computation time of the simulation, only the domain of A-E-O-H-Á-É-Ó-H’ is used to create the numerical simulation model.

In this study, since the climate data and the heat fluxes, presented in Figure 10, are only available for 1 year [33], then the climate data and the heat fluxes are repeated 50 times to generate a 50-year data. The undisturbed temperature of the ground is considered to be 4.41 °C, equal to the annual average of the equivalent temperature at the ground surface. The diameter of borehole is set to be 0.112 m. The ground material is considered to be homogenous without any groundwater. The thermal properties of the ground including thermal conductivity, specific heat capacity and density are considered to be $3\mathrm{W}/(\mathrm{m}·\mathrm{K})$, $750\mathrm{J}/(\mathrm{kg}·\mathrm{K})$ and $2500\mathrm{kg}/{\mathrm{m}}^3$ [49]. In this study, the latent heat associated with the phase change does not taken into account, so the ground temperature below 0 °C only indicates possible ground freezing.

Furthermore, the surface boundary conditions of the ground; i.e., the area of A-B-C-D, are simplified by an equivalent temperature, $T_{eq} (°\mathrm{C})$, and an equivalent heat transfer coefficient, $h_{eq} (\mathrm{W}/({\mathrm{m}}^2·\mathrm{K}))$. If ${\tilde{h}}_c (\mathrm{W}/({\mathrm{m}}^2·\mathrm{K}))$ is the annual mean value of the convective heat transfer coefficient and ${\tilde{h}}_r (\mathrm{W}/({\mathrm{m}}^2·\mathrm{K}))$ is the annual mean value of radiation heat transfer coefficient, the values of $T_{eq}$ and $h_{eq}$ are calculated as:

$$h_{eq}={\tilde{h}}_c+{\tilde{h}}_r$$

(30)

$$T_{eq}=\frac{{\tilde{h}}_c·T_{ambient}+{\tilde{h}}_r·T_{sky}+\alpha ·I}{h_{eq}}$$

(31)

From the climate data for Östersund [33], the value of ${\tilde{h}}_c$ is $13.7\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$ with the standard deviation of 4.8 $\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$ and the value of ${\tilde{h}}_r$ is $4.3\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$ with the standard deviation of 0.4 $\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$. For details about calculating ${\tilde{h}}_c$ and ${\tilde{h}}_r$, the reader is referred to [29]. Furthermore, the boundary conditions at bottom and side surfaces of the ground volume of A-B-C-D-Á-B’-Ć-D’ are set to be adiabatic.

The results of the average temperature evolution at the borehole walls over 50-year period associated with different inlet temperatures of the fluid in the HHP system and three different ${\mathrm{COP}}_{\mathrm{h}}$ are presented in

Table 4. The results related to the spatial average temperature evolution at the borehole walls over 50-year period for different inlet temperatures of the HHP system and three different ${\mathrm{COP}}_{\mathrm{h}}$s.

${\mathit{T}}_{\mathit{i}\mathit{n}}$ (°C)	${\mathit{T}}_{\mathit{l}\mathit{i}\mathit{m}\mathit{i}\mathit{t}}$ (°C)	COPh = 3				COPh = 4				COPh = 5
		${\overline{\mathit{q}}}_{\mathit{s}-\mathit{y}\mathit{e}\mathit{a}\mathit{r}}$ $(\mathrm{W}/\mathrm{m})$	${\mathit{T}}_{\mathit{a}\mathit{v}\mathit{e}}$ (°C)	${\mathit{T}}_{\mathit{a}\mathit{m}\mathit{p}}$ (°C)	${\mathit{T}}_{\mathit{c}\mathit{h}\mathit{a}\mathit{n}\mathit{g}\mathit{e}}$ (°C)	${\overline{\mathit{q}}}_{\mathit{s}-\mathit{y}\mathit{e}\mathit{a}\mathit{r}}$ $(\mathrm{W}/\mathrm{m})$	${\mathit{T}}_{\mathit{a}\mathit{v}\mathit{e}}$ (°C)	${\mathit{T}}_{\mathit{a}\mathit{m}\mathit{p}}$ (°C)	${\mathit{T}}_{\mathit{c}\mathit{h}\mathit{a}\mathit{n}\mathit{g}\mathit{e}}$ (°C)	${\overline{\mathit{q}}}_{\mathit{s}-\mathit{y}\mathit{e}\mathit{a}\mathit{r}}$ $(\mathrm{W}/\mathrm{m})$	${\mathit{T}}_{\mathit{a}\mathit{v}\mathit{e}}$ (°C)	${\mathit{T}}_{\mathit{a}\mathit{m}\mathit{p}}$ (°C)	${\mathit{T}}_{\mathit{c}\mathit{h}\mathit{a}\mathit{n}\mathit{g}\mathit{e}}$ (°C)
4	6	1.5	4.9	2.1	0.5	1.5	4.9	2.1	0.5	1.4	4.9	2.1	0.5
6	6	1.2	4.8	2.0	0.4	1.1	4.8	2.1	0.4	1.1	4.8	2.1	0.4
8	8	0.9	4.7	1.9	0.3	0.8	4.7	2.0	0.3	0.8	4.7	2.0	0.3
10	8	0.6	4.6	1.9	0.2	0.6	4.6	1.9	0.2	0.5	4.6	2.0	0.2
12	10	0.4	4.5	1.8	0.1	0.3	4.5	1.9	0.1	0.3	4.5	1.9	0.1

. The results are: the inlet temperature of the fluid, $T_{in} (°\mathrm{C})$, the temperature for turning on the harvesting operation, $T_{limit} (°\mathrm{C})$, the annual average heat rate injected/extracted to or from a single borehole, ${\overline{q}}_{s-year} (\mathrm{W}/\mathrm{m})$, the annual spatial average temperature at the borehole walls during 50 years, $T_{ave} (°\mathrm{C})$, the annual amplitude temperature at the borehole walls, $T_{amp} (°\mathrm{C})=(T_{max} (°\mathrm{C})-T_{min} (°\mathrm{C}))/2$, and the annual average temperature change at the borehole wall from the first year to 50th year, $T_{change} (°\mathrm{C})$.

Figure 12a shows hourly BTES load for one year which is calculated based on heat fluxes of the HHP system, see Figure 10, and Equation (27). Moreover, Figure 12b shows the hourly variation of the spatial average temperature evolution at the borehole walls over 50-year period associated with $T_{in}=6 °\mathrm{C}$ and $T_{limit}=6 °\mathrm{C}$ and ${\mathrm{COP}}_{\mathrm{h}}=4$.

As can be seen from Table 4, the annual average temperature change at the borehole walls over time is above 0 °C. Furthermore, as it can be seen from Figure 12b, the temperature evolution over time follows an increasing trend. The temperature increase is significant for the first two years and then slowly trends towards zero with time. The temperature increase at the borehole walls indicates that the annual heat injected into the ground during harvesting period is higher than the sum of the annual heat extracted from the ground during anti-icing period and the heat losses from the BTES to the surrounding ground. This means that BTES with 20 borehole and 200 m depth for $T_{in}$ varies between $4°\mathrm{C}$ and $12°\mathrm{C}$ and ${\mathrm{COP}}_{\mathrm{h}}$ between 3 and 5 can operate reliably over 50 years of operation with a slight improvement compared to the first year.

To understand the thermal process of the temperature change at the borehole walls, the main model of the BTES is separated into four sub-models using the superposition principle, see Figure 13.

Figure 13a is related to the sub-model which is dealing with the initial temperature of the BTES. The initial temperature is equal to the annual average value of the equivalent temperature, see Equation (31). The remaining thermal process in the ground can be divided into three parts originating from: (i) an annual average value, $q_{average} (\mathrm{W}/{\mathrm{m}}^2)$, (ii) a periodical variation, $q_{periodical} (\mathrm{W}/{\mathrm{m}}^2)$, due to the yearly cycling of the injected/extracted heat from the borehole round the average value and (iii) the deviation initial temperature, $T_{dev} (\mathrm{K})$. In Figure 13c, it is assumed that the periodic process exists from time $-\infty$ to $+\infty$. To get the initial temperature ${\overline{T}}_{eq}$ in the whole ground by superimposing all sub-models in Figure 13, the initial temperature at time t = 0 for the last sub-model, $T_{dev} (\mathrm{K})$, should be equal to the ground temperature from the periodical sub-model, Figure 13c, but with negative sign. Furthermore, In Figure 13d, the influence of the variations in the outdoor temperature during the year is neglected. This variation will only influence the ground temperature a few meters down from the surface. $q_{average}$ is the main part of the heat fluxes which leads to the temperature change at the borehole walls. The value of $q_{average}$ can be obtained using Equation (28). $q_{periodical}$ gives a zero contribution to the annual average borehole temperature and $T_{dev}$ only leads to a small disturbance during the first few years. However, in this study, this temperature disturbance is assumed to be negligible. Hence, the results of sub-model (b) can be used as the representative of the annual average temperature changes at the borehole walls. Comparing the results related to the annual average temperature changes at the borehole walls over 50-year period obtained from sub-model (b) and the results presented in Table 4 shows that maximum relative difference between the results is 3%.

6. Conclusions

By assuming the initial geometry of the HHP system in Table 3, the BTES geometry from Figure 11 and the climate data of Östersund, it was found from this study that:

• An increase in the inlet fluid temperature resulted in a decrease in the number of hours of the slippery conditions. Increasing the inlet fluid temperature from 4 °C to 20 °C led to approximately a 65% decrease in the number of hours of the slippery conditions on the road surface.
• Considering a constant air temperature as a criterion for turning on/off the harvesting operation in the HHP system, an increase in the inlet fluid temperature resulted in a decrease in the harvested solar energy. However, by considering a constant inlet fluid temperature, the harvested solar energy depended on the air temperature which was used for turning on/off the harvesting process.
• A decrease in the inlet fluid temperature and an increase in the air temperature as a criterion for turning on/off the harvesting operation resulted in an increase in the temperature reduction on the road surface during harvesting period. Regardless of when the harvesting operation start, the inlet fluid temperature of 4 °C resulted in a more than 5 °C temperature reduction on the road surface during harvesting period.
• The BTES, consisted of 20 boreholes and 200 m depth, led to a reasonable hourly temperature fluctuation at the borehole walls between $2.7°\mathrm{C}$ to $7.1°\mathrm{C}$. Analyzing the temperature increase/decrease at the borehole walls showed that the BTES can operate consistently and reliably in the long term for all cases that the annual harvested solar energy in summer is higher than the annual required energy for anti-icing the road surface in winter. For the inlet fluid temperature lower than or equal to $12°\mathrm{C}$, the harvested solar energy was higher than the required energy for anti-icing the road surface.

This study was a basic concept of using the HHP system and BTES for harvesting solar energy and anti-icing the road surface. Coupling the HHP system to the BTES including heat pumps as well as the details related to the design of BTES are needed to be investigated in further studies. In addition, it is interesting to optimize the number of boreholes using the thermal loads during the harvesting and anti-icing the road surfaces and evaluate the seasonal energy storage using the exergy performance. Furthermore, it is needed to investigate the harvesting solar energy and anti-icing operations with varying inlet fluid temperatures and the air velocity passing on the surface of the road. In this study, it was assumed that the temperature variation does not affect the thermal properties of the materials and the fluid. Further studies are needed to simulate the HHP and the BTES with varying thermal properties, corresponding to the temperature of the materials and the fluid. The main goal of constructing the HHP system is to use solar energy to provide a sustainable solution to mitigate the slippery conditions. However, the sustainability should be investigated by considering wider aspects such as the initial investment and operation costs, life span of the road with HHP system versus conventional roads, reducing the traffic accidents and increasing the road accessibility, and finally, the environmental impacts of the HHP system.