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

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 kWh/(^m2·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 kWh/(^m2·year) with a standard deviation of 7 kWh/(^m2·year) and the mean value of the required energy is 119 kWh/(^m2·year) with a standard deviation of 73 kWh/(^m2·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×¹⁰⁻⁴ 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 (kg/^m2) is the mass balance of the water on the road surface, ^m˙″ _con (kg/(^m2·s)) is the condensation rate and ^m˙″ _evp (kg/(^m2·s)) is the evaporation rate, the mass balance of water on the road surface is obtained as:

d_mmoisture″dt=_m˙con″−_m˙evp″

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, _qsnow (W/m²), related to the heat capacity of snow and temperature change is considered in the heat balance:

1. conductive heat from ground and pipes:

_qcond=−λ·∇T

2. convective heat flow from the ambient air

_qconv=_hc·(_Tambient−_Tsurface)

3. sensible heat from rain

_qrain=_m˙rain·_cp−water·(_Tambient−_Tsurface)

4. sensible heat from snow

_qsnow=_m˙snow·_cp−snow·(_Tambient−_Tsurface)

5. long-wave radiation

_qlw=ε·σ·(_Tsky4−_Tsurface4)

6. short-wave radiation

_qsw=α·I

7. latent heat of evaporation and condensation

_qevp/con=_he·β·(_vambient−_vs)

where λ (W/(m·K)) is the thermal conductivity of the road materials, T (K) is the temperature, _Tambient (K) is the ambient air temperature, _Tsurface (K) is the road surface temperature, _hc (W/(^m2·K)) is the convective heat transfer coefficient, ε (-) is the emissivity of the surface, σ (W/(^m2·^K4)) is the Stefan-Boltzmann constant, _Tsky (K) is the sky temperature, α (-) is the solar absorptivity of the surface, I (W/^m2) is the solar irradiation, _he (J/kg) is the latent heat of evaporation of water, β (m/s) is the moisture transfer coefficient, _vambient (kg/^m3) is the humidity by the volume of the ambient air, _vs (kg/^m3) is the humidity by the volume of the saturated air at the surface temperature and m˙ (kg/(^m2·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:

_qcond+_qconv+_qrain+_qsnow+_qlw+_qsw+_qevp/con=0

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 _Req−pipe (m·K/W) . _Req−pipe consists of three serially coupled thermal resistance, namely: (i) the surface heat transfer resistance at inner pipe surface, _RPWS (m·K/W) , (ii) the pipe resistance, _Rp (m·K/W) , and (iii) the thermal resistance between the outer pipe wall and an additional annulus, _Ri.j (m·K/W) . The value of _RPWS is calculated using Reynolds, Prandtl and Nusselt dimensionless numbers, see Equations (10) to (13).

Re=2·_ρf·_vf·_rinnerμf

Pr=_μf·_cp,fλf

Nu={4 if Re≤23000.023·^Re0.8·^Pr1/3 if Re>2300

_RPWS=1π·_λf·Nu

where Re (-) is the Reynolds number, Pr (-) is the Prandtl number, Nu (-) is the Nusselt number, _ρf (kg/^m3) is the density of fluid, _vf (m/s) is the average velocity of fluid through a cross section of pipe, _rinner (m) is the inner radius of the pipe, _uf (kg/(m·s)) is the dynamic viscosity of fluid, _cp,f (J/(kg·K)) is the specific heat capacity of fluid and _λf (W/(m·K)) is the thermal conductivity of fluid. Furthermore, _Rp is calculated as:

_Rp=ln(_routerrinner)2·π·_λpipe

where _rinner (m) and _router (m) are respectively the inner and outer radii of the pipe and _λpipe (W/(m·K)) is the thermal conductivity of the pipe material. Finally, _Ri,j is calculated as:

_Ri,j=ln(_ri,jrouter)2·π·_λi,j

where _ri,j (m) is the radius of the additional annulus surrounding the pipe and _λi,j (W/(m·K)) is the thermal conductivity of the surrounding materials. It should be noted that _ri,j>_router . The total thermal resistance between the fluid and the embedding materials, _Req−pipe , is obtained from Equation (16). Moreover, the equivalent temperature surrounding the pipe, _Teq−pipe (K) , corresponds to the average ambient temperature of the embedding materials close to the pipe segment is calculated as Equation (17).

_Req−pipe=_Rp+_Ri,j+_RPWS

_Teq−pipe=_Tf−_qi,j·_Req−pipe

where _Tf (K) is the temperature of fluid, circulating through the pipe and _qi,j (W/^m2) 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).

_vf·π·_rinner ²·_ρf·_cp,f·∂_Tf(x)∂x+_Tf(x)−_Teq−pipe(x)_Req−pipe=0

where x (m) is the longitudinal length coordinate along the fluid motion, _Tf(x) and _Teq−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:

1. The transversal heat transfer flow for each 2D section of the HHP system is calculated independently of the adjacent section.

2. The fluid temperature distribution is calculated using a local quasi steady state condition. The condition is valid for each time step.

3. The fluid temperature distribution is calculated using a step-wise sequence solution for all 2D sections of the HHP systems.

4. 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 _L1 , _L2 , … _Ln .

For the pipe section between the labels of n and n + 1, the section length is _Ln , the inlet temperature of the fluid at time t is _Tf,nt and the outlet temperature is _Tf,n+1t . Considering the steady state longitudinal fluid temperature distribution along the pipe, the analytical solution to obtain the outlet fluid temperature is as:

_Tf,n+1t=_{Teq−pipe,nt}+(_Tf,nt−_{Teq−pipe,nt})·^{e−(_Ln/_ln)}

where _ln (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 _ln is calculated as:

_ln=_Req−pipe·_vf·π·_rinner ²·_ρf·_cf

In this study, the inlet temperature of fluid, _Tf,1t (K), is a known parameter. By applying Equations (19) and (20), it is possible to calculate the _Tf,2t (K). The value of _Tf,2t 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 W/^m2 . 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.

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:

1. Case 1: the ambient air temperature was 25 °C and the initial temperature was 20 °C.

2. 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. 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 °C±0.2 °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 °C and 25 °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. 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, _Tsurface (K), is below the dew-point temperature, _Tdew (K), the water will be formed on the road surface due to the condensation. In addition, if _Tsurface is below the freezing temperature of water, _Tfreezing (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 _Tdew , when _Tsurface<_Tfreezing . The mentioned criteria for running the heating system could be written as:

{_Tsurface<_TdewTsurface<_Tfreezing

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, _Er (kWh/(m²·year)), is calculated as:

_Er=1c·L·_{∫t=01 year}(_Tf,outt−_Tf,int)·_vf·π·_rinner ²·_ρf·_cp,f·dt

where c (m) is the distance between pipes, L (m) is the pipe length, _Tf,int (K) is the inlet temperature of the fluid and _Tf,outt (K) is the outlet temperature of the fluid, _vf (m/s) is the velocity of the fluid, _rinner (m) is the inner radius of the pipe, _ρf (kg/^m3) is the density of the fluid and _cp,f (J/(kg·K)) is the specific heat capacity of the fluid.

Furthermore, the number of hours of the slippery condition on the road surface, _tslippery (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:

_tslippery=_{∫01 year}f·dt (f=1 if (_Tsurface<_TdewTsurface<_Tfreezing) else f=0)

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 °C [26]. Considering a lower air temperature such as 4 °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 °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 °C and for about 10% of this period the air temperature is above 16 °C [33]. In this study, five different air temperatures of 6 °C , 8 °C , 10 °C , 12 °C and 16 °C are arbitrary selected as a criterion to start harvesting solar energy. If _Tlimit (°C) is one of the mentioned five temperatures, then the criterion to start harvesting solar energy can be written as:

_Tair>_Tlimit

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

_Tdecrease=_Tunheated−_Tharv

where _Tunheated (°C) is the surface temperature of an unheated road and _Tharv (°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, _Eh (kWh/(^m2·year)) , can be calculated as:

_Eh=_{∫t=01 year}(_Tout−r−_Tin−r)·_vf·π·_rinner ²·_ρf·_cp,f·dt_cr·_Lr (if _Tair>_Tlimit)

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 _Tsurface higher than (but close to) _Tdew when _Tsurface<_Tfreezing , see Section 3.4. Considering that the maximum air temperature in Östersund is 25 °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, _Tin (°C) , from 4 °C to 20 °C leads to an increase in the annual required energy for anti-icing the road surface, _Er , from 81.2 kWh/(^m2·year) to 96.3 kWh/(^m2·year) . Higher _Tin causes higher surface temperature and consequently results in more energy consumption. However, an increase in the _Tin from 4 °C to 16 °C results in a decrease in the annual harvested solar energy, _Eh , on average from 271.7 kWh/(^m2·year) to 34.7 kWh/(^m2·year) . Interesting point is that for _Tin=20 °C , the HHP system is mainly working as a heater rather than a solar collector, so as 44.4 kWh/(^m2·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 _Tlimit , see Equation (24). By increasing _Tlimit from 6 °C to 16 °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 _Tin and a lower value of _Tlimit can result in heating the road surface during harvesting period, rather than harvesting solar energy. For example, considering _Tin=20 °C and _Tlimit=6 °C results in consuming 140 kWh/(^m2·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, _Er=96.3 kWh/(^m2·year) , when _Tin=20 °C . To harvest maximum possible solar energy from the road surface, it is essential to use a suitable combination of _Tin and _Tlimit . The combinations of 4 °C≤_Tin≤6 °C and _Tlimit=6 °C , 8 °C≤_Tin≤10 °C and _Tlimit=8 °C , _Tin=12 °C and _Tlimit=10 °C , _Tin=16 °C and _Tlimit=12 °C as well as _Tin=20 °C and _Tlimit=16 °C lead to the maximum values of _Eh associated with each inlet fluid temperature, see Figure 9a. The maximum harvested solar energy during summer is 352.1 kWh/(^m2·year) which is related to _Tin=4 °C and _Tlimit=6 °C . It is worth mentioning that for 4 °C≤_Tin≤12 °C , the harvested solar energy can be higher than the required energy for anti-icing the road surface. However, for _Tin≥16 °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, _Tout , versus _Tin . During both harvesting and anti-icing periods, an increase in _Tin will results in an increase in _Tout . However, the absolute difference between the inlet and outlet temperatures, |_Tin−_Tout| , during harvesting and heating periods is different. During harvesting period and for a constant value of _Tlimit , an increase in _Tin results in harvesting less solar energy, which in turn leads to a lower value of |_Tin−_Tout| . On the contrary, during anti-icing period, an increase in _Tin results in consuming more energy for heating the road surface, which in turn leads to a higher value of |_Tin−_Tout| .

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

Although an increase in _Tin causes a decrease in _Tdecrease , 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, _tslippery . From Figure 9d, by increasing _Tin from 4 °C to 20 °C, the value of _tslippery reduces from 284 h to 93 h. It should be noted that the decrease of _tslippery follows a convex curve; i.e., as the value of _Tin increases the number of slippery hours approximates to the preceding value of _tslippery . It is worth mentioning that using Equation (1) and considering the density of ice as 917 kg/^m3 [31], the maximum height of remaining ice on the road surface decreases from 1.48 mm to 0.98 mm as the value of _Tin 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 _Tin≤12 °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, _qp(W/^m2) , is extracted from the hybrid 3D numerical simulation model in Section 4. The heat flux of _qp 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, _qp=0 W/^m2 . 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 _qp 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, _{q¯p_year} (W/^m2) is the annual average of heat fluxes at the road surface. As can be seen, the value of _{q¯p_year} is positive and varying from 30.9 W/^m2 for _Tin=4 °C and _Tlimit=6 °C to 5.2 W/^m2 for _Tin=12 °C and _Tlimit=10 °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 _COPh is the coefficient of performance of the heat pump during anti-icing period, the heat flux adjustment factor for anti-icing, _Hf , can be defined as [44]:

_Hf=_COPh−1_COPh

By multiplying the required heat flux of the HHP system by _Hf , the extracted heat from the BTES can be estimated during the anti-icing process. In this study, for simplicity of the simulation, three constant _COPh s of 3, 4 and 5 are used [45]. It is worth mentioning that using a constant _COPh may affect some errors in the hourly simulation results, so by a 1 °C decrease in the entering fluid temperature to the heat pump, the value of _COPh 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, _qs(W/m) , can be obtained using the harvested/anti-icing heat fluxes as:

_qs=_Wr·_Lrdtotal·_qp×f ({f=1 during harvestingf=_Hf during anti−icing)

where _Lr (m) is the length of the HHP system, _Wr (m) is the width of the HHP system and _dtotal (m) is the total required depth of borehole. In this study, _Lr=50 m and _Wr=3.5 m . The value of _dtotal can be determined using the following equation as [47]:

_dtotal=_qh·_Rb+_qy·_R10y+_qm·_R1m+_qh·_R6hTm−(_Tg+_Tp)(_dtotal=N×_dbore)

where N (-) is the number of boreholes, _dbore (m) is the depth of each borehole, _qh (W) is the maximum hourly ground peak load, _qm (W) is the average monthly ground load during the month when _qh occurs, _qy (W) is the net heat load of the ground and _Rb (m·K/W) is the effective borehole thermal resistance. _R10y (m·K/W) , _R1m (m·K/W) and _R6h (m·K/W) represent, respectively, the effective ground thermal resistance of yearly, monthly, and hourly ground load components. _Tm (°C) is the mean fluid temperature in the borehole, _Tg (°C) is the undisturbed ground temperature and _Tp (°C) is the temperature penalty. For more details of the parameters, the reader is referred to [47,48].

Considering _Rb=0.1 m·K/W and using the average heat loads related to the different fluid temperatures from 4 °C to 12 °C and three different _COPh s of 3, 4 and 5, the value of _dtotal is calculated as 3898 m. In this simulation, the depth of each borehole, _dbore , 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×_dbore . 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 W/(m·K) , 750 J/(kg·K) and 2500 kg/^m3 [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, _Teq (°C) , and an equivalent heat transfer coefficient, _heq (W/(^m2·K)) . If _h˜c (W/(^m2·K)) is the annual mean value of the convective heat transfer coefficient and _h˜r (W/(^m2·K)) is the annual mean value of radiation heat transfer coefficient, the values of _Teq and _heq are calculated as:

_heq=_h˜c+_h˜r

_Teq=_h˜c·_Tambient+_h˜r·_Tsky+α·I_heq

From the climate data for Östersund [33], the value of _h˜c is 13.7 W/(^m2·K) with the standard deviation of 4.8 W/(^m2·K) and the value of _h˜r is 4.3 W/(^m2·K) with the standard deviation of 0.4 W/(^m2·K) . For details about calculating _h˜c and _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 _COPh are presented in Table 4. The results are: the inlet temperature of the fluid, _Tin (°C) , the temperature for turning on the harvesting operation, _Tlimit (°C) , the annual average heat rate injected/extracted to or from a single borehole, _q¯s−year (W/m) , the annual spatial average temperature at the borehole walls during 50 years, _Tave (°C) , the annual amplitude temperature at the borehole walls, _Tamp (°C)=(_Tmax (°C)−_Tmin (°C))/2 , and the annual average temperature change at the borehole wall from the first year to 50th year, _Tchange (°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 _Tin=6 °C and _Tlimit=6 °C and _COPh=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 _Tin varies between 4 °C and 12 °C and _COPh 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, _qaverage (W/^m2) , (ii) a periodical variation, _qperiodical (W/^m2) , due to the yearly cycling of the injected/extracted heat from the borehole round the average value and (iii) the deviation initial temperature, _Tdev (K) . In Figure 13c, it is assumed that the periodic process exists from time −∞ to +∞ . To get the initial temperature _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, _Tdev (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. _qaverage is the main part of the heat fluxes which leads to the temperature change at the borehole walls. The value of _qaverage can be obtained using Equation (28). _qperiodical gives a zero contribution to the annual average borehole temperature and _Tdev 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:

1. 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.

2. 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.

3. 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.

4. The BTES, consisted of 20 boreholes and 200 m depth, led to a reasonable hourly temperature fluctuation at the borehole walls between 2.7 °C to 7.1 °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 °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.

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]."]

[Image omitted. See PDF.] and the hybrid 3D numerical simulation model developed in this study for two different cases (a) the temperature of asphalt pavement at the depth of 2.5 cm and (b) the outlet temperature of fluid."]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

Material	Thermal Conductivity (W/(m·K))	Density (kg/m3)	Specific Heat Capacity (J/(kg·K))
concrete pavement	1.62	2339	896
Pipe	0.42	932	2300
Fluid (water)	0.61	1000	4181

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

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·ri,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

Tin (°C)	Tlimit (°C)	COPh = 3				COPh = 4				COPh = 5
		q¯s−year (W/m)	Tave (°C)	Tamp (°C)	Tchange (°C)	q¯s−year (W/m)	Tave (°C)	Tamp (°C)	Tchange (°C)	q¯s−year (W/m)	Tave (°C)	Tamp (°C)	Tchange (°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

Author Contributions

R.M. is a PhD student, C.-E.H. is the main supervisor and P.J. is the co-supervisor.

Funding

This work was funded by the Norwegian Public Road Administration and Chalmers University of Technology.

Acknowledgments

The first author would like to acknowledge Mohamad Kharseh and Josef Johnsson from Chalmers University of Technology for their interest in the work and fruitful discussions.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

Symbols

c Distance between pipes (m)

_cp specific heat capacity (J/(kg·K))

_COPh Coefficient of performance (-)

_dbore Depth of borehole (m)

D Depth of the embedded pipe (m)

DIA Diameter (m)

E Energy (kWh/(^m2·year)

_hc Convective heat transfer coefficient (W/^m2·K)

_he Latent heat of water evaporation (kJ/kg)

_heq Equivalent heat transfer coefficient (W/^m2·K)

_hr Radiation heat transfer coefficient (W/^m2·K)

I Solar irradiation (W/^m2)

l Characteristic length (m)

L Length (m)

m Mass per square meter ( kg/^m2 )

m˙ Mass rate per square meter (kg/(^m2·s))

N Number of boreholes (-)

Nu Nusselt number (-)

Pr Prandtl number (-)

q Heat flux (W/^m2)

_q¯p-year Annual average of heat flux of road (W/^m2)

_qs Heat rate of borehole (W/m)

_q¯s-year Annual average of heat rate of borehole (W/m)

r Radius of the pipe (m)

R Thermal resistance (m·K/W)

_Rs Surface thermal resistance (^m2·K/W)

Re Reynolds number (-)

t Time (s)

T Temperature (K)

_T¯eq Annual average of equivalent temperature (K)

_Tlimit Criterion for turning on harvesting operation (K)

_Tmax Annual max. temperature at borehole wall (K)

_Tmin Annual min. temperature at borehole wall (K)

v Velocity (m/s)

_V˙f Fluid flow rate (L/min)

W Width (m)

x Longitudinal length coordinate (m)

Greek symbols

α Solar absorptivity (-)

β Surface moisture transfer coefficient (m/s)

ε Emissivity coefficient (-)

ρ Density (kg/^m3)

σ Stephan-Boltzmann constant (-)

λ Thermal conductivity (W/(m·K))

v Humidity by the volume (kg/^m3)

u Dynamic viscosity (kg/(m·s))

Subscripts

cond Conductive heat

conv Convective heat

dev Deviation

dew Dew-point

eq Equivalent temperature/resistance

evp/con Evaporation and condensation

f Fluid

g Ground

h Harvested energy

i,j Additional annulus surrounding the pipe

in Inner/Inlet

lw Long-wave radiation

m Month/Mean

n Number of sections

out Outer/Outlet

p Pipe/Penalty

PWS Inner pipe wall surface

r Required energy, road

s Surface, Storage

sw Short-wave radiation

y year

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

References

1. Andersson, A.; Chapman, L. The use of a temporal analogue to predict future traffic accidents and winter road conditions in Sweden. Meteorol. Appl. 2011, 18, 125-136.

2. Strandroth, J.; Rizzi, M.; Olai, M.; Lie, A.; Tingvall, C. The effects of studded tires on fatal crashes with passenger cars and the benefits of electronic stability control (ESC) in Swedish winter driving. Accid. Anal. Prev. 2012, 45, 50-60.

3. Arvidsson, A.K. The Winter Model-A new way to calculate socio-economic costs depending on winter maintenance strategy. Cold Reg. Sci. Technol. 2017, 136, 30-36.

4. Wåhlin, J.; Klein-Paste, A. A salty safety solution. Phys. World 2017, 30, 27-30.

5. Asensio, E.; Ferreira, V.J.; Gil, G.; García-Armingol, T.; López-Sabirón, A.M.; Ferreira, G. Accumulation of de-icing salt and leaching in Spanish soils surrounding roadways. Int. J. Environ. Res. Public Health 2017, 14, 1498.

6. Anuar, L.; Amrin, A.; Mohammad, R.; Ourdjini, A. Vehicle accelerated corrosion test procedures for automotive in Malaysia. MATEC Web Conf. 2017, 90, 01040.

7. Knudsen, F.; Natanaelsson, K.; Arvidsson, A.; Kärki, O.; Jacobsen, Á.; Guðmundsson, G.F.; Nonstad, B.; Reitan, K.M. Vintertjeneste i de Nordiske land. Statusrapport 2014; NVF-Rapporter: Norge, Norway, 2014. (In Norwegian)

8. Liu, T.; Pan, Q.; Sanchez, J.; Sun, S.; Wang, N.; Yu, H. Prototype Decision Support System for Black Ice Detection and Road Closure Control. IEEE Intell. Transp. Syst. Mag. 2017, 9.

9. Gáspár, L.; Bencze, Z. Salting Route Optimization in Hungary. Transp. Res. Procedia 2016, 14, 2421-2430.

10. Elvik, R. Does the influence of risk factors on accident occurrence change over time? Accid. Anal. Prev. 2016, 91, 91-102.

11. Shaopeng, W.; Mingyu, C.; Jizhe, Z. Laboratory Investigation into Thermal Response of Asphalt Pavements as Solar Collector by Application of Small-Scale Slabs. Appl. Therm. Eng. 2011, 31, 1582-1587.

12. Cheng, Y.C.; Tao, J.L.; Jiao, Y.B.; Guo, Q.L.; Li, C. Influence of Diatomite and Mineral Powder on Thermal Oxidative Ageing Properties of Asphalt. Adv. Mater. Sci. Eng. 2015, 2015, 947834.

13. Wang, X.; Gu, X.; Ni, F.; Deng, H.; Dong, Q. Rutting resistance of porous asphalt mixture under coupled conditions of high temperature and rainfall. Constr. Build. Mater. 2018, 174, 293-301.

14. Mirzanamadi, R. Ice Free Roads Using Hydronic Heating Pavement with Low Temperature: Thermal Properties of Asphalt Concretes and Numerical. Licentiate Thesis, Degree of Licentiate of Engineering, Chalmers University of Technology, Gothenburg, Sweden, 2017.

15. Wang, H.; Jasim, A.; Chen, X. Energy harvesting technologies in roadway and bridge for different applications-A comprehensive review. Appl. Energy 2018, 212, 1083-1094.

16. Ceylan, H.; Gopalakrishnan, K.; Kim, S.; Cord, W. Heated Transportation Infrastructure Systems: Existing and Emerging Technologies. In Proceedings of the 12th International Symposium on Concrete Roads, Prague, Czech Republic, 23-26 September 2014.

17. Sundberg, J. Halkfria Vägar-Studier av Typfall Solvärme och Värmelagring för Miljöanpassad Halkbekämpning; Trafikverket: Götenburg, Sweden, 2014. (In Swedish)

18. Eugster, W.J. Road and Bridge Heating Using Geothermal Energy. In Overview and Examples. In Proceedings of the European Geothermal Congress, Unterhaching, Germany, 30 May-1 June 2007.

19. Pahud, D. Serso, Stockage Saisonnier Solaire pour le Dégivrage d'un Pont. Rapport Final; Office fédéral de l'énergie: Bern, Switzerland, 2007. (In French)

20. Pahud, D. BRIDGESIM: Simulation Tool for the System Design of Bridge Heating for Ice Prevention with Solar Heat Stored in a Seasonal Ground Duct Store, User Manual; SUPSI: Lugano, Switzerland, 2008.

21. Abbasi, M. Non-Skid Winter Road, Investigation of Deicing System by Considering Different Road Profiles. Master's Thesis, Chalmers University of Technology, Gothenburg, Sweden, 2013.

22. Gao, L.; Zhao, J.; Tang, Z. A Review on Borehole Seasonal Solar Thermal Energy Storage. Energy Procedia 2015, 70, 209-218.

23. Johnsson, J. Winter Road Maintenance Using Renewable Thermal. Licentiate Thesis, Degree of Licentiate of Engineering, Chalmers University of Technology, Gothenburg, Sweden, 2017.

24. Liu, X. Development and Experimental Validation of Simulation of Hydronic Snow Melting Systems for Bridges. Ph.D. Thesis, Oklahoma State University, Stillwater, OK, USA, 2005.

25. Javed, S. Design of Ground Source Heat Pump Systems-Thermal Modelling and Evaluation of. Licentiate Thesis, Degree of Licentiate of Engineering, Chalmers University of Technology, Gothenburg, Sweden, 2010.

26. Mirzanamadi, R.; Hagentoft, C.-E.; Johansson, P.; Johnsson, J. Anti-icing of road surfaces using Hydronic Heating Pavement with low temperature. Cold Reg. Sci. Technol. 2018, 145, 106-118.

27. Li, B.; Wu, S.P.; Xiao, Y.; Pan, P. Investigation of heat-collecting properties of asphalt pavement as solar collector by a three-dimensional unsteady model. Mater. Res. Innov. 2015, 19, S1-172-S1-176.

28. Mirzanamadi, R. Utilizing Solar Energy for Anti-Icing Road Surfaces Using Hydronic Heating Pavement with Low Temperature. Ph.D. Thesis, Chalmers University of Technology, Gothenburg, Sweden, 2019.

29. Mirzanamadi, R.; Hagentoft, C.-E.; Johansson, P. An analysis of hydronic heating pavement to optimize the required energy for anti-icing. Appl. Therm. Eng. 2018, 144, 278-290.

30. Mirzanamadi, R.; Hagentoft, C.-E.; Johansson, P. Coupling a Hydronic Heating Pavement to a Horizontal Ground Heat Exchanger for Harvesting Solar Energy and Heating Road Surfaces. Appl. Energy 2018, submit.

31. Hagentoft, C.-E. Introduction to Building Physics; 1:7; Studentlitteratur: Lund, Sweden, 2001; ISBN 9789144018966.

32. Karlsson, H. Embedded water-based surface heating part 1: Hybrid 3D numerical model. J. Build. Phys. 2010, 33, 357-391.

33. Meteotest. Meteonorm: Meteonorm, Global Meteorological Database. Handbook Part II: Theory, Version 6.1; Meteotest: Bern, Switzerland, 2010.

34. Mirzanamadi, R.; Johansson, P.; Grammatikos, S.A. Thermal properties of asphalt concrete: A numerical and experimental study. Constr. Build. Mater. 2018, 158, 774-785.

35. Andersland, O.B.; Ladanyi, B. Frozen Ground Engineering, 2nd ed.; John Wiley & Sons: Hoboken, NJ, USA, 2004; ISBN 978-0-471-61549-1.

36. Bergman, T.L.; Incropera, F.P.; Lavine, A.S.; Dewitt, D.P. Introduction to Heat Transfer, 6th ed.; Wiley and Sons: New York, NY, USA, 2011.

37. Eppelbaum, L.; Kutasov, I.; Pilchin, A. Applied Geothermics. In Applied Geothermics; Springer: Berlin/Heidelberg, Germany, 2014; ISBN 978-3-642-34022-2.

38. Robertson, E.C. Report of Thermal Properties of Rocks; United State Department of the Interior Geological Survey: Reston, VA, USA, 1988.

39. Bohne, D.; Fischer, S.; Obermeier, E. Thermal Conductivity, Density, Viscosity, and Prandtl-Numbers of Ethylene Glycol-Water Mixtures. Ber. Bunsenges. Phys. Chem. 1984, 88, 739-742.

40. Pan, P.; Wu, S.P.; Xiao, Y.; Liu, G. A Review on Hydronic Asphalt Pavement for Energy Harvesting and Snow Melting. Renew. Sustain. Energy Rev. 2015, 48, 624-634.

41. Bobes-Jesus, V.; Pascual-Muñoz, P.; Castro-Fresno, D.; Rodriguez-Hernandez, J. Asphalt Solar Collectors: A Literature Review. Appl. Energy 2013, 102, 962-970.

42. Bayer, P.; de Paly, M.; Beck, M. Strategic optimization of borehole heat exchanger field for seasonal geothermal heating and cooling. Appl. Energy 2014, 136, 445-453.

43. You, T.; Wu, W.; Shi, W.; Wang, B.; Li, X. An overview of the problems and solutions of soil thermal imbalance of ground-coupled heat pumps in cold regions. Appl. Energy 2016, 177, 515-536.

44. Law, Y.L.E.; Dworkin, S.B. Characterization of the effects of borehole configuration and interference with long term ground temperature modelling of ground source heat pumps. Appl. Energy 2016, 179, 1032-1047.

45. Banks, D. An Introduction to Thermalgeology: Ground Source Heating and Cooling, 2nd ed.; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2012; ISBN 978-1-4051-7061-1.

46. Spitler, J.D. Glhepro-A Design Tool for Commercial Building Ground Loop Heat Exchangers. In Proceedings of the Fourth International Heat Pumps in Cold Climates Conference, Aylmer, QC, Canada, 17-18 August 2000.

47. Monzó, P.; Bernier, M.; Acuña, J.; Mogensen, P. A monthly based bore field sizing methodology with applications to optimum borehole spacing. ASHRAE Trans. 2016, 122, 111-126.

48. Bernier, M.A.; Chahla, A. Long-Term Ground-Temperature Changes in Geo-Exchange Systems. ASHRAE Trans. 2008, 114, 342-350.

49. Javed, S. Thermal Modelling and Evaluation of Borehole Heat Transfer. Ph.D. Thesis, Chalmers University of Technology, Gothenburg, Sweden, 2012.