Energy Analysis of a Dual-Source Heat Pump Coupled with Phase Change Materials

: Installation costs of ground heat exchangers (GHEs) make the technology based on ground-coupled heat pumps (GCHPs) less competitive than air source heat pumps for space heating and cooling in mild climates. A smart solution is the dual source heat pump (DSHP) which switches between the air and ground to reduce frosting issues and save the system against extreme temperatures a ﬀ ecting air-mode. This work analyses the coupling of DSHP with a ﬂat-panel (FP) horizontal GHE (HGHE) and a mixture of sand and phase change materials (PCMs). From numerical simulations and considering the energy demand of a real building in Northern Italy, di ﬀ erent combinations of heat pumps (HPs) and trench backﬁll material were compared. The results show that PCMs always improve the performance of the systems, allowing a further reduction of the size of the geothermal facility. Annual average heat ﬂux at FP is four times higher when coupled with the DSHP system, due to the lower exploitation. Furthermore, the enhanced dual systems are able to perform well during extreme weather conditions for which a sole air source heat pump (ASHP) system would be unable either to work or perform e ﬃ ciently. Thus, the DSHP and HGHE with PCMs are robust and resilient alternatives for air conditioning.


Introduction
Ground-coupled heat pumps (GCHPs) are among the general renewable technologies used for space heating, cooling and hot water supply for buildings, especially in cold climates. Significant reductions in energy consumption, 30-70% in heating mode and 20-50% in cooling mode, have been reported when compared with conventional air-conditioning systems [1]. The ground heat exchanger (GHE), installed vertically or horizontally to exchange heat with the soil, is one of the main components of the GCHP system. Although GCHPs are more efficient than conventional air-to-water heat pumps [2], the investment cost of the GHE is high and make GCHPs less competitive than the more widespread air source heat pumps (ASHPs), especially in mild climates. Thus, despite having great potential in the future for many countries [3], the presence of GCHPs in countries with mild climates is currently not very widespread. Some novel and shallow horizontal GHE (HGHE) solutions take advantage of a better heat transfer from advanced shapes, achieving an energy performance similar to that obtained by the more thermally stable, though more expensive, vertical configurations. In this regard, new developments of very shallow geothermal systems combined with different natural and cheap backfilling materials are currently under assessment at several sites around Europe [4]. However, it seems HGHEs hold some drawbacks regarding land use [5] and payback still remains too long to justify the initial investment [6]. Yet, given the existing potential and resources, a major were paraffin A8 (PCM1) and paraffin A24 (PCM2) [22], respectively. According to their manufacturer [22], all these organic PCMs are comprised of long chain molecules, usually with a carbon backbone. The longer the chain, the higher the melting point. They were selected according to the meteorological conditions, specifically those related to minimum and maximum temperatures, in the site where heating and cooling demand were set, and ground temperatures trends of preliminary simulations. Their main thermal properties are shown in Table 1, as well as those proposed elsewhere [23] for wet sand (trench) and wet clay (surrounding soil). Volume ratios of PCM1 and PCM2 and sand into the trench were calculated considering a sand porosity reference (40%), and according to the gross building energy demand for heating (2/3) and cooling (1/3) throughout the year for the location studied.

Building Energy Demand and Heat Pump Efficiency
The energy requirements of the building of TekneHub Laboratory (University of Ferrara) for a one-year period were considered as the heating and cooling demand of the system. TekneHub is located in the city of Ferrara (44°50′ N, 11°37′ E) in Northern Italy. Taking this building as a reference case, EnergyPlus building energy simulation program [24] was used to estimate thermal energy loads  The flat shape of the flat-panel solution allows an easier coupling with PCMs than piping solutions of other HGHEs, such as straight pipes, slinky coils or baskets. The panel was 1.5 m high, providing a total heat transfer surface of 3 m 2 per metre trench length. The trench, which was 2.8 m deep and 0.6 m wide, was filled with a mixture of wet sand and two different PCMs, whilst the surrounding soil was assumed as wet clay. Considering cooling and heating modes, the PCMs used were paraffin A8 (PCM1) and paraffin A24 (PCM2) [22], respectively. According to their manufacturer [22], all these organic PCMs are comprised of long chain molecules, usually with a carbon backbone. The longer the chain, the higher the melting point. They were selected according to the meteorological conditions, specifically those related to minimum and maximum temperatures, in the site where heating and cooling demand were set, and ground temperatures trends of preliminary simulations. Their main thermal properties are shown in Table 1, as well as those proposed elsewhere [23] for wet sand (trench) and wet clay (surrounding soil). Volume ratios of PCM1 and PCM2 and sand into the trench were calculated considering a sand porosity reference (40%), and according to the gross building energy demand for heating (2/3) and cooling (1/3) throughout the year for the location studied.

Building Energy Demand and Heat Pump Efficiency
The energy requirements of the building of TekneHub Laboratory (University of Ferrara) for a one-year period were considered as the heating and cooling demand of the system. TekneHub is located in the city of Ferrara (44 • 50 N, 11 • 37 E) in Northern Italy. Taking this building as a reference Energies 2020, 13, 2933 4 of 17 case, EnergyPlus building energy simulation program [24] was used to estimate thermal energy loads (q t , W/m 3 ) every hour for the year 2015 ( Figure 2), by using local weather data collected from a meteorological station operating at TekneHub. The building envelope is well isolated according to recent Italian regulations on the energy performance of buildings. The air conditioning system consists of two air-to-air rooftop heat pumps with a capacity of 40 kW each. Further details of the building components and load calculation can be found in a study by Bottarelli and colleagues [25].  Figure 2), by using local weather data collected from a meteorological station operating at TekneHub. The building envelope is well isolated according to recent Italian regulations on the energy performance of buildings. The air conditioning system consists of two air-to-air rooftop heat pumps with a capacity of 40 kW each. Further details of the building components and load calculation can be found in a study by Bottarelli and colleagues [25]. Regarding the heat pump, heating and cooling efficiencies were measured by the coefficient of performance (COP) and the energy efficiency ratio (EER), respectively. These ratios basically depend on the temperature difference between the condenser and evaporator, which is strictly related to the heat source temperatures for fixed indoor temperature set points (20 °C in winter, 26 °C in summer). The performance curves for the heat pump used by Nam et al. in their research [26] have been used here: (1) For simplicity, evaporator temperature (Te) and condenser temperature (Tc) are assumed to be equal to the temperature of the sources (air, ground heat exchanger) in the following.

Numerical Modelling and Boundary Conditions
The main aim here is to calculate the variation in the temperature distribution of the ground source, due to the extraction carried out by a flat-panel HGHE in coupling with a GCHP or a DSHP, and with or without a PCM mixture filling the trench. This kind of problem is governed by the equation of heat conduction in a solid (Equation (3)), and it has been numerically solved by using the commercial software COMSOL Multiphysics [27], in which some equations were implemented to characterise the PCM mixture according to the equivalent heat capacity method, as further detailed.
A two-dimensional (2D) computational domain which represents a cross section of the shallow GHE installation is considered, assuming the following: uniform ground temperature distribution along the z-axis, a small temperature change between the inlet and outlet along the exchanger, a negligible thermal resistance of the flat-panel, and no thermal stratification of the working fluid. Thus, the average temperature at the flat-panel is considered representative of the working fluid temperature. The simplification of the wall thermal resistance of the GHE could compromise a direct comparison with the air heat exchanger. However, it should be considered that the flat-panel prototype is actually made with walls 1.6 mm thick in high-density polyethylene (HDPE) with thermal conductivity around 0.48 W/mK). Therefore, its wall thermal resistance is around 0.0033 K/W, which is comparable with the fouling thermal resistance that affects an air heat exchanger, and Regarding the heat pump, heating and cooling efficiencies were measured by the coefficient of performance (COP) and the energy efficiency ratio (EER), respectively. These ratios basically depend on the temperature difference between the condenser and evaporator, which is strictly related to the heat source temperatures for fixed indoor temperature set points (20 • C in winter, 26 • C in summer). The performance curves for the heat pump used by Nam et al. in their research [26] have been used here: COP = 0.002T 2 e + 0.0597T e + 2.9028 (1) For simplicity, evaporator temperature (T e ) and condenser temperature (T c ) are assumed to be equal to the temperature of the sources (air, ground heat exchanger) in the following.

Numerical Modelling and Boundary Conditions
The main aim here is to calculate the variation in the temperature distribution of the ground source, due to the extraction carried out by a flat-panel HGHE in coupling with a GCHP or a DSHP, and with or without a PCM mixture filling the trench. This kind of problem is governed by the equation of heat conduction in a solid (Equation (3)), and it has been numerically solved by using the commercial software COMSOL Multiphysics [27], in which some equations were implemented to characterise the PCM mixture according to the equivalent heat capacity method, as further detailed.
A two-dimensional (2D) computational domain which represents a cross section of the shallow GHE installation is considered, assuming the following: uniform ground temperature distribution along the z-axis, a small temperature change between the inlet and outlet along the exchanger, a negligible thermal resistance of the flat-panel, and no thermal stratification of the working fluid. Thus, the average temperature at the flat-panel is considered representative of the working fluid temperature. The simplification of the wall thermal resistance of the GHE could compromise a direct comparison with the air heat exchanger. However, it should be considered that the flat-panel prototype is actually made with walls 1.6 mm thick in high-density polyethylene (HDPE) with thermal Energies 2020, 13, 2933 5 of 17 conductivity around 0.48 W/mK). Therefore, its wall thermal resistance is around 0.0033 K/W, which is comparable with the fouling thermal resistance that affects an air heat exchanger, and even an order lower if recent techno-polymers (with thermal conductivity greater than 2 W/mK) are used instead of a standard HDPE. Figure 3 depicts the domain and its boundary conditions. Due to symmetry conditions, only half of the domain is simulated. Figure 4 provides ground surface temperature at TekneHub Laboratory for 2015, which is highly correlated with air temperature. Maximum, mean and minimum temperatures values are 34 • C, 16 • C and 2 • C, respectively. This time series (T s ) has been set on the top of the domain (first-type boundary condition) instead of developing a more sophisticated energy balance equation at the ground surface (third-type boundary condition), since both methodologies lead to similar ground and GHE wall temperature values beyond 0.8 m deep, as reported by Bortoloni and co-workers [5]. A constant temperature of 16.7 • C was set at the bottom of the domain, corresponding to the annual average temperature measured (T b ).
Energies 2020, 13, x FOR PEER REVIEW 5 of 17 even an order lower if recent techno-polymers (with thermal conductivity greater than 2 W/mK) are used instead of a standard HDPE. Figure 3 depicts the domain and its boundary conditions. Due to symmetry conditions, only half of the domain is simulated. Figure 4 provides ground surface temperature at TekneHub Laboratory for 2015, which is highly correlated with air temperature. Maximum, mean and minimum temperatures values are 34 °C, 16 °C and 2 °C, respectively. This time series (Ts) has been set on the top of the domain (first-type boundary condition) instead of developing a more sophisticated energy balance equation at the ground surface (third-type boundary condition), since both methodologies lead to similar ground and GHE wall temperature values beyond 0.8 m deep, as reported by Bortoloni and co-workers [5]. A constant temperature of 16.7 °C was set at the bottom of the domain, corresponding to the annual average temperature measured (Tb).  The flat-panel GHE is simplified as a boundary heat source/sink (second-type boundary condition), for which hourly heat flux (qFP, W/m 2 ) is calculated from Equation (4): where qt is the thermal load per unit volume of building space (W/m 3 ) at time t and described in Section 2.2; SFP is the heat transfer surface of a flat-panel HGHE per unit of length of the trench (2 Energies 2020, 13, x FOR PEER REVIEW 5 of 17 even an order lower if recent techno-polymers (with thermal conductivity greater than 2 W/mK) are used instead of a standard HDPE. Figure 3 depicts the domain and its boundary conditions. Due to symmetry conditions, only half of the domain is simulated. Figure 4 provides ground surface temperature at TekneHub Laboratory for 2015, which is highly correlated with air temperature. Maximum, mean and minimum temperatures values are 34 °C, 16 °C and 2 °C, respectively. This time series (Ts) has been set on the top of the domain (first-type boundary condition) instead of developing a more sophisticated energy balance equation at the ground surface (third-type boundary condition), since both methodologies lead to similar ground and GHE wall temperature values beyond 0.8 m deep, as reported by Bortoloni and co-workers [5]. A constant temperature of 16.7 °C was set at the bottom of the domain, corresponding to the annual average temperature measured (Tb).  The flat-panel GHE is simplified as a boundary heat source/sink (second-type boundary condition), for which hourly heat flux (qFP, W/m 2 ) is calculated from Equation (4): where qt is the thermal load per unit volume of building space (W/m 3 ) at time t and described in Section 2.2; SFP is the heat transfer surface of a flat-panel HGHE per unit of length of the trench (2  The flat-panel GHE is simplified as a boundary heat source/sink (second-type boundary condition), for which hourly heat flux (q FP , W/m 2 ) is calculated from Equation (4): where q t is the thermal load per unit volume of building space (W/m 3 ) at time t and described in Section 2.2; S FP is the heat transfer surface of a flat-panel HGHE per unit of length of the trench (2 sides times 1.5 m 2 /m); l f is the load factor calculated as the ratio between the building space volume (m 3 ) and the length of the trench (m). The parameter l f provides a smart way to parametrically link the required length of the flat-panel with building energy demand (volume-to-length ratio, [21]). From the parametric study developed by Bottarelli and co-workers [21], it was concluded that the DSHP solution offered the best combination of efficiency and opportunity to reduce installation cost when l f is between 10 m 3 /m and 15 m 3 /m. In the present study, l f values varied within the range 10 m 3 /m to 30 m 3 /m, with 10 m 3 /m and 20 m 3 /m as reference values for the GCHP and DSHP, respectively. Finally, an adiabatic condition was set to the rest of the side boundaries of the domain.
The problem was solved numerically by using COMSOL Multiphysics [27]. Using the finite element method, a preliminary grid independence study was carried out. The case with the presence of PCMs and the dual source functionality was considered to evaluate the average temperature at the flat-panel in an unsteady state for a period of 40 days and six different meshes. Newton's iterative method was selected to solve the system of equations at each time-step. Absolute and relative tolerances were set to 10 −5 • C and 10 −3 , respectively. Figure 5 depicts all resulting values; relative changes in the estimation of average flat-panel temperature were smaller than 1%. A final mesh consisting of around 7000 triangular elements (Figure 3b), finer around the GHE trench, was finally selected to simulate all cases.
Energies 2020, 13, x FOR PEER REVIEW 6 of 17 sides times 1.5 m 2 /m); lf is the load factor calculated as the ratio between the building space volume (m 3 ) and the length of the trench (m). The parameter lf provides a smart way to parametrically link the required length of the flat-panel with building energy demand (volume-to-length ratio, [21]). From the parametric study developed by Bottarelli and co-workers [21], it was concluded that the DSHP solution offered the best combination of efficiency and opportunity to reduce installation cost when lf is between 10 m 3 /m and 15 m 3 /m. In the present study, lf values varied within the range 10 m 3 /m to 30 m 3 /m, with 10 m 3 /m and 20 m 3 /m as reference values for the GCHP and DSHP, respectively. Finally, an adiabatic condition was set to the rest of the side boundaries of the domain.
The problem was solved numerically by using COMSOL Multiphysics [27]. Using the finite element method, a preliminary grid independence study was carried out. The case with the presence of PCMs and the dual source functionality was considered to evaluate the average temperature at the flat-panel in an unsteady state for a period of 40 days and six different meshes. Newton's iterative method was selected to solve the system of equations at each time-step. Absolute and relative tolerances were set to 10 −5 °C and 10 −3 , respectively.  When the mixture of sand and PCMs is considered within the trench, heat transfer is simplified to heat conduction within an equivalent solid mixture, according to the equivalent heat capacity method. Thus, the equivalent mixture density and thermal conductivity are calculated as volume weighted averages according to component volume proportions, ri (Equations (5) and (6)). In contrast with other studies [16,21], the specific heat of the mixture was calculated as a mass weighted averaged (Equation (7)). These mixture properties are then calculated as follows: When the mixture of sand and PCMs is considered within the trench, heat transfer is simplified to heat conduction within an equivalent solid mixture, according to the equivalent heat capacity method. Thus, the equivalent mixture density and thermal conductivity are calculated as volume weighted averages according to component volume proportions, r i (Equations (5) and (6)). In contrast with other studies [16,21], the specific heat of the mixture was calculated as a mass weighted averaged (Equation (7)). These mixture properties are then calculated as follows: Energies 2020, 13, 2933 7 of 17 where r i and Y i are the volume and mass fractions of each i component, respectively. These ratios are related as follows: In the present study, two PCMs are considered, such that n = 2 (i = 0 for sand, 1 for PCM1 and 2 for PCM2). H i (T) and D i (T) are step and normalised Dirac's pulse functions, respectively, used to represent the evolution of the different physical properties during the phase change, similarly to the approach reported by Bottarelli et al. [16]. Specifically, c P i (T) includes the specific heat capacity in each phase (c S P i , c L P i ) and latent heat fusion (h i ) for PCM's i component. The distribution of latent heat is approximated by using the D i (T) function, in such a way that the total heat per unit volume released during the phase transformation equals the latent heat. A transitional interval of 8 • C centred at the melting temperature was considered. As sand porosity is not greater than 40%, the sum of the PCMs' volume ratios (r 1 + r 2 ) is limited and set to 0.40; only in a single case was this value exceeded. In addition, according to the ratio of the periods of heating and cooling demand (Figure 2), r 1 and r 2 were set equal to 0.25 and 0.15, respectively.
Furthermore, a virtual control of the boundary condition at flat-panel was programmed in COMSOL in order to simulate the switching between air and ground sources in the DSHP case, as described by Bottarelli and colleagues [21]. The flowchart in Figure 6 describes the performance routine. In winter, the ground achieves better working conditions when its temperature is higher than outdoor air temperature. Ground source is used in the model when outdoor air temperature is smaller than 5 • C (T air limit ), and also when ground source temperature (T g ) exceeds air temperature by an onset temperature differential, ∆T. According to the guidelines for the maximum benefits of DSHPs reported by Bottarelli and colleagues [21], the corresponding temperature differential was around 5 • C for the range of variation of l f (from 10 to 30 in this study). Both rules (T air limit = 5 • C; ∆T = 5 • C) aim to preserve the ground source for the late winter, because the exploitation of the high ground temperature at the winter beginning would lead to quick depletion of the underground thermal energy storage.
In the present study, two PCMs are considered, such that n = 2 (i = 0 for sand, 1 for PCM1 and 2 for PCM2). Hi(T) and Di(T) are step and normalised Dirac's pulse functions, respectively, used to represent the evolution of the different physical properties during the phase change, similarly to the approach reported by Bottarelli et al. [16]. Specifically, includes the specific heat capacity in each phase ( , ) and latent heat fusion (hi) for PCM's i component. The distribution of latent heat is approximated by using the Di(T) function, in such a way that the total heat per unit volume released during the phase transformation equals the latent heat. A transitional interval of 8 °C centred at the melting temperature was considered. As sand porosity is not greater than 40%, the sum of the PCMs' volume ratios (r1 + r2) is limited and set to 0.40; only in a single case was this value exceeded. In addition, according to the ratio of the periods of heating and cooling demand (Figure 2), r1 and r2 were set equal to 0.25 and 0.15, respectively.
Furthermore, a virtual control of the boundary condition at flat-panel was programmed in COMSOL in order to simulate the switching between air and ground sources in the DSHP case, as described by Bottarelli and colleagues [21]. The flowchart in Figure 6 describes the performance routine. In winter, the ground achieves better working conditions when its temperature is higher than outdoor air temperature. Ground source is used in the model when outdoor air temperature is smaller than 5 °C (Tair limit), and also when ground source temperature (Tg) exceeds air temperature by an onset temperature differential, ΔT. According to the guidelines for the maximum benefits of DSHPs reported by Bottarelli and colleagues [21], the corresponding temperature differential was around 5 °C for the range of variation of lf (from 10 to 30 in this study). Both rules (Tair limit = 5 °C; ΔT = 5 °C) aim to preserve the ground source for the late winter, because the exploitation of the high ground temperature at the winter beginning would lead to quick depletion of the underground thermal energy storage.
An operating threshold is finally set for the ground source temperature (Tg lim = −2 °C) to prevent working fluid in the flat-panel from freezing. In summer, ground source is on when the difference between air and ground source temperatures is higher than the onset differential.
Air source temperature was defined as that of outdoor air (year 2015) at an hourly time scale. The temperature of the ground source was estimated by the model every time step and reported at an hourly scale as the average temperature at the flat-panel boundary.   An operating threshold is finally set for the ground source temperature (T g lim = −2 • C) to prevent working fluid in the flat-panel from freezing. In summer, ground source is on when the difference between air and ground source temperatures is higher than the onset differential.
Air source temperature was defined as that of outdoor air (year 2015) at an hourly time scale. The temperature of the ground source was estimated by the model every time step and reported at an hourly scale as the average temperature at the flat-panel boundary.
The initial temperature distribution in the solid domain was derived from experimental measurements of ground temperature at different depths ( Figure 7). However, all models were initially run for a full year to achieve an initial condition already affected by their specific yearly exploitation.
Energies 2020, 13, x FOR PEER REVIEW 8 of 17 The initial temperature distribution in the solid domain was derived from experimental measurements of ground temperature at different depths ( Figure 7). However, all models were initially run for a full year to achieve an initial condition already affected by their specific yearly exploitation.

Case Studies
Different case studies have been numerically simulated in unsteady heat transfer conditions in order to quantify the potential benefits of using both the dual system, specifically a DSHP coupled with the flat-panel (FP) HGHE, and a mixture of sand and PCM as backfill material. Water to water GCHP with sand as backfill material was taken as a reference case. Table 2 describes all the cases considered. Load factor (lf) was increased in relation to the reference case (DP00), for which lf equals 10 m 3 /m (building cubic metre over metre FP length), in order to assess the energy performance and capability of the enhanced systems by using smaller heat transfer surfaces of the flat-panel and PCMs. Indeed, numerical simulations previously carried out showed that higher values of lf for the DP00 leaded to temperatures several degrees below zero, which are not compliant with common operating conditions of this kind of facility. Furthermore, the effect of a reduction of trench width (wt) was also analysed in case DP11#, as well as higher volume PCM ratios (r1 = 0.50 and r2 = 0.20) and higher thermal PCM conductivity (k1 = k2 = 1.2 W/mK) in case DP11+.
In all cases, the novel flat-panel HGHE is coupled with the GCHP or DSHP. The bias between the experiments and the numerical model implemented in COMSOL of this HGHE is less than 1 °C on average, as reported by Bottarelli and colleagues [21]. When the DSHP is on, the value of the onset temperature differential (ΔT) was set to 5 °C.

Results and Discussion
All models were run for a two-year period by repeating weather data of the year 2015 and under the same boundary conditions, in order to set the typical initial conditions of this kind of facility and

Case Studies
Different case studies have been numerically simulated in unsteady heat transfer conditions in order to quantify the potential benefits of using both the dual system, specifically a DSHP coupled with the flat-panel (FP) HGHE, and a mixture of sand and PCM as backfill material. Water to water GCHP with sand as backfill material was taken as a reference case. Table 2 describes all the cases considered. Load factor (l f ) was increased in relation to the reference case (DP00), for which l f equals 10 m 3 /m (building cubic metre over metre FP length), in order to assess the energy performance and capability of the enhanced systems by using smaller heat transfer surfaces of the flat-panel and PCMs. Indeed, numerical simulations previously carried out showed that higher values of l f for the DP00 leaded to temperatures several degrees below zero, which are not compliant with common operating conditions of this kind of facility. Furthermore, the effect of a reduction of trench width (w t ) was also analysed in case DP11#, as well as higher volume PCM ratios (r 1 = 0.50 and r 2 = 0.20) and higher thermal PCM conductivity (k 1 = k 2 = 1.2 W/mK) in case DP11+. In all cases, the novel flat-panel HGHE is coupled with the GCHP or DSHP. The bias between the experiments and the numerical model implemented in COMSOL of this HGHE is less than 1 • C on Energies 2020, 13, 2933 9 of 17 average, as reported by Bottarelli and colleagues [21]. When the DSHP is on, the value of the onset temperature differential (∆T) was set to 5 • C.

Results and Discussion
All models were run for a two-year period by repeating weather data of the year 2015 and under the same boundary conditions, in order to set the typical initial conditions of this kind of facility and to ensure that all models reproduced their stationary trend according to their different exploitations. Thus, Figure 8a shows flat-panel temperatures for the second simulation year cases DP00, DP01, DP10 and DP11, together with outdoor air and undisturbed ground temperature (at a depth of 1.7 m). Similarly, Figure 8b shows cases DP11, DP11*, DP11# and DP11+. As it can be inferred from the series of undisturbed ground temperature, despite ground thermal exploitation and unlike deep geothermal systems, thermal drift is avoided by using this shallow HGHE.
°C, such that anti-freeze usage is needed. As a consequence, dual strategy leaves a warmer ground at the beginning of summer which is less advantageous for this period. On the contrary, the lower temperature of GCHP systems makes the condition disadvantageous in wintertime, but more profitable in summer. When PCMs are used in the backfill material (DP01 and DP11) a similar behaviour is found, although the performance of the system is clearly lower in wintertime (minimum temperatures of −2.4 °C and 1.7 °C). These results are mainly due to the low thermal conductivity of the mixture in comparison with that of the sand. However, despite this drawback, the performance is improved in summertime, according to the higher cold energy stored in the ground, as maximum temperatures of 24.9 °C and 25.7 °C are achieved for DP01 and DP11, respectively. In Figure 8b, DP11* and DP11# minimum temperatures are −0.3 °C and -0.6 °C in wintertime, respectively, whilst 1.3 °C is that of DP11+. In summertime, the maximum temperatures are 26.1 °C, 26.2 °C and 24.7 °C for DP11*, DP11# and DP11+, respectively.
Overall, dual systems are able to perform well during extreme weather conditions (very low or very high outdoor air temperatures) for which a sole ASHP system would be unable either to work or perform efficiently. This makes the dual system a robust alternative in Southern European countries in which weather conditions are expected to become more severe in the future, with higher inter-annual increase in summer temperatures and low variability in current winter temperatures [28]. Details of the annual trend are shown at the beginning of the year in Figure 9, and by the middle of the year in Figure 10; that is, in winter and summer, respectively. When comparing, it should be noted that case DP11* has a higher lf factor (30 m 3 /m), DP11# also includes a narrow trench and consequently a smaller quantity of PCMs, whilst DP11+ has a standard trench (60 cm wide) and load factor (20 m 3 /m), but uses PCMs with higher thermal conductivity.
In Figure 9a, the steady decrease of flat-panel temperature in January due to a continuous ground exploitation carried out by the GCHP (DP00 and DP01) contrasts with the nearly constant trend shown by the DSHP (DP10, DP11, DP11*). DP11+ shows the promptest celerity in recovering Flat-panel temperature for DSHP cases (DP10 and DP11) is higher than for GCHP cases (DP00 and DP01), since the dual system is able to switch to the air when ground temperature is lower in winter (Figure 8a). In this way, the lowest flat-panel temperatures for DP10 and DP11 cases are 2.1 • C and 1.7 • C, respectively. This issue may avoid the use of anti-freeze (e.g., propylene glycol) in the secondary loop of the system and reduce or exclude all frosting problems at the air heat exchanger. In contrast, the lowest flat-panel temperatures for GCHP cases DP00 and DP01 are −1.3 • C and −2.4 • C, such that anti-freeze usage is needed. As a consequence, dual strategy leaves a warmer ground at the beginning of summer which is less advantageous for this period. On the contrary, the lower temperature of GCHP systems makes the condition disadvantageous in wintertime, but more profitable in summer. When PCMs are used in the backfill material (DP01 and DP11) a similar behaviour is found, although the performance of the system is clearly lower in wintertime (minimum temperatures of −2.4 • C and 1.7 • C). These results are mainly due to the low thermal conductivity of the mixture in comparison with that of the sand. However, despite this drawback, the performance is improved in summertime, according to the higher cold energy stored in the ground, as maximum temperatures of 24.9 • C and 25.7 • C are achieved for DP01 and DP11, respectively. In Figure 8b, DP11* and DP11# minimum temperatures are −0.3 • C and −0.6 • C in wintertime, respectively, whilst 1.3 • C is that of DP11+. In summertime, the maximum temperatures are 26.1 • C, 26.2 • C and 24.7 • C for DP11*, DP11# and DP11+, respectively.
Overall, dual systems are able to perform well during extreme weather conditions (very low or very high outdoor air temperatures) for which a sole ASHP system would be unable either to work or perform efficiently. This makes the dual system a robust alternative in Southern European countries in which weather conditions are expected to become more severe in the future, with higher inter-annual increase in summer temperatures and low variability in current winter temperatures [28].
Details of the annual trend are shown at the beginning of the year in Figure 9, and by the middle of the year in Figure 10; that is, in winter and summer, respectively. When comparing, it should be noted that case DP11* has a higher l f factor (30 m 3 /m), DP11# also includes a narrow trench and consequently a smaller quantity of PCMs, whilst DP11+ has a standard trench (60 cm wide) and load factor (20 m 3 /m), but uses PCMs with higher thermal conductivity.    In Figure 9a, the steady decrease of flat-panel temperature in January due to a continuous ground exploitation carried out by the GCHP (DP00 and DP01) contrasts with the nearly constant trend shown by the DSHP (DP10, DP11, DP11*). DP11+ shows the promptest celerity in recovering the exploitation in comparison with all other DSHPs (Figure 9b). Cases with PCMs show better performance than those without (Figure 10a), whilst the aforementioned behaviour of DP11+ is confirmed also in summertime (Figure 10b). Therefore, this last property seems to be the key factor in using PCMs: not to penalise the improvement of energy storage with a lower thermal conductivity of the backfill material.
In Figure 11, the resulting performance of the GCHP (DP00, DP01) and DSHP (DP10, DP11) are depicted in terms of COP (wintertime, Figure 11a) and EER (summertime, Figure 11b), according to the previous Equations (1) and (2). GCHP cases show better COP values than DSHP cases at the beginning of winter, when the ground temperature is very high due to the previous heating up in summer. However, the reverse happens late in winter, when the ground has been deeply exploited by the GCHP operation, and partially saved by the DSHP mode. This inversion does not happen in summertime (EER), since the lowest ground temperature performs better during the whole summer. depicted in terms of COP (wintertime, Figure 11a) and EER (summertime, Figure 11b), according to the previous Equations (1) and (2). GCHP cases show better COP values than DSHP cases at the beginning of winter, when the ground temperature is very high due to the previous heating up in summer. However, the reverse happens late in winter, when the ground has been deeply exploited by the GCHP operation, and partially saved by the DSHP mode. This inversion does not happen in summertime (EER), since the lowest ground temperature performs better during the whole summer.  A summary of the energy exchange and performance of the different systems in winter, summer and the whole year is shown in Tables 3-5, respectively. The values of thermal (Qt) and electrical (Qe) energy exchanged per unit length of trench (kWh/m) are given too. In wintertime (Table 3), DP10 and DP11 perform similarly, and better than DP00 and DP01. Specifically, DP11 shows a higher COP than DP01 (3.44 against 3.24), and also DP10 (COP value of 3.43) compared to DP00 (COP value of 3.26). Furthermore, DP11* and DP11# achieve very similar performance if compared with DP11, but with a significant reduction of the HGHE size (l f of 30m 3 /m against 20 m 3 /m), and a narrower trench for case DP11# (40 cm against 60 cm).
Due to the much shorter use of the ground heat exchanger carried out by the DSHP cases (DP10, DP11, DP11*, DP11# and DP11+), very high values of heat flux per metre of trench are obtained, that range from more than three times for case DP11+ up to almost five times for case DP11*, when compared with GCHP. This is mainly related to the shorter ground exploitation of the DSHP cases, for which heat transfer takes from 500 h (DP11#) up to 720 h (DP11+).
This seems to reflect a good performance of the dual system even if coupled with a shorter HGHE, and that PCMs can improve the system performance only if their thermal conductivities are higher. Therefore, the common low thermal conductivity of PCMs attenuates their potential good performance in wintertime related to their underground thermal energy storage.
In comparison with the cheapest ASHP which shows an overall COP value around 3.31, the best value of a DSHP, around 3.46 (DP11#) does not seem to justify the investment of PCMs and GHE. However, it should be noted that the air temperature drops below 3 • C for 747 h over 5228 h of heating time, that is, 14.2% of the working period, in which frosting issues are very common at this location. As a consequence, the ASHP performance given here should be considered overestimated up to a value around 15%, as reported by Dongellini and colleagues [29].
In summertime, DP01 shows the best performance (EER value of 4.37), as justified by the ground cooling occurring for long thermal exploitation in wintertime, and the presence of PCMs. The best DSHP case is DP11# with an EER value of 4.02, resulting from the 4.24 and 3.80 values in ground and air exploitation, respectively. It is worth noting that DP11# uses a GHE which is three times shorter than the one used by case DP00. The ASHP performance is the lowest (EER value of 3.61 value) and a more interesting gap is evident, especially related to a larger exploitation of the GHE. In summer season, the overall needs of air conditioning cover a period of 1304 h; when the DSHP is used, the ground exploitation goes from a minimum of 488 h (37.4%, DP11#) to a maximum of 635 h (48.7%, DP11+). The average heat flux occurring at the flat-panel GHE (q FP_ave ) shows a minimum of 50.4 W per length metre of trench in GCHP cases (DP00, DP01) and a maximum of 142.8 W/m in DSHP case (DP11#), which demonstrates the high performance of the flat-panel. q FP_ave values are quite comparable in winter and summer for DSHP cases, while winter values halved those in summer for GCHP cases. This difference is always attributable to the exploitation time of the ground, unchanged for DSHP cases (500-700 h) unlike GCHP cases (with more than 5000 h in wintertime against 1300 h in summertime).
As expected from the annual trend of flat-panel temperature, all DSHP cases show higher efficiencies than GCHP cases. Yet, they still perform better than the ASHP. Regarding the heat flux per metre of trench, DSHP cases provide values which are several times those given by GCHP cases. During the whole year, DP00 and DP01 show COP/EER values of 3.42, the ASHP of 3.37 and the DSHP from 3.48 to 3.56. Overall, the average heat flux occurring at the flat-panel, coupled with a DSHP system, is three to five times higher than that of a GCHP.
Finally, an analysis about how the PCM behaves in the different cases was carried out. Thus, Figure 12 shows the solid phase time series in terms of mass fraction for the trench. In winter, a more efficient use of PCM1 is made by DP01 (Figure 12a) than by DP11 (Figure 12b). Hence, PCM1 solidifies reaching a peak of 100% (DP01), while it is of 80% in DP11. A slight improvement is achieved by DP11*, while DP11+ gets slightly lower results (Figure 12c,d). In summertime, the fraction of PCM2 that becomes liquid is very similar in all cases. This result could also suggest the advisability of using backfill materials in the trench with higher thermal conductivity. To complete the results depiction, a sequence of the thermal field at 41.333 days of simulation time is presented in Figure 13 for different cases. The first two images (DP00, DP01) show the large exploitation carried out by the GCHP in comparison with the DSHP; no relevant differences are shown when PCMs are used in coupling with a GCHP. More interesting differences are detectable among DSHP cases. In the sequence DP10, DP11 and DP11*, the ground temperature rises in the domain from the bottom to the top, whilst the flat-panel temperature decreases. Temperature distribution for DP11 and DP11# cases are quite similar, while DP11+ shows the highest temperature values among all. To complete the results depiction, a sequence of the thermal field at 41.333 days of simulation time is presented in Figure 13 for different cases. The first two images (DP00, DP01) show the large exploitation carried out by the GCHP in comparison with the DSHP; no relevant differences are shown when PCMs are used in coupling with a GCHP. More interesting differences are detectable among DSHP cases. In the sequence DP10, DP11 and DP11*, the ground temperature rises in the domain from the bottom to the top, whilst the flat-panel temperature decreases. Temperature distribution for DP11 and DP11# cases are quite similar, while DP11+ shows the highest temperature values among all. exploitation carried out by the GCHP in comparison with the DSHP; no relevant differences are shown when PCMs are used in coupling with a GCHP. More interesting differences are detectable among DSHP cases. In the sequence DP10, DP11 and DP11*, the ground temperature rises in the domain from the bottom to the top, whilst the flat-panel temperature decreases. Temperature distribution for DP11 and DP11# cases are quite similar, while DP11+ shows the highest temperature values among all.

Conclusions
In this study, the energy performance of the coupling of a DSHP and a novel flat-panel HGHE with a mixture of sand and PCM backfill material into the trench was analysed numerically. To the authors' knowledge, this combination, used both for space heating and cooling, has not been considered or studied previously in the scientific literature. A methodology based on energy numerical simulations, in which different source heat pumps (GCHP, ASHP and DSHP) and configurations (with and without PCMs, and different GHEs and trench sizes) were compared under the same real operating conditions, was developed to check the suitability of the system. The main conclusions of the study are as follows.
Numerical simulations show that the overall performance of a ground source heat pump is higher, although not very significantly, than the one given by an air source heat pump. The use of shallow HGHEs avoids thermal drift of ground temperature and shows significant values of heat flux value per unit length of the exchanger, even when compared to vertical boreholes. The flat-panel shows an annual average heat flux of around 33 W/m when it is coupled with a GCHP, and more than 140 W/m with a DSHP. The dual-source functionality optimises the performance of a heat pump and drastically allows the reduction of the HGHE length to up to several times that of a full GCHP, causing a significant reduction in installation costs. Moreover, the value estimated for the lowest working fluid temperature precludes the use of glycol, with benefits in terms of costs and environmental impact. The use of PCMs as backfilling material can improve the overall energy performance only if a significant increase of their thermal conductivity is carried out, as the common low value of the pure paraffins' conductivity reduces the benefit of their high latent heat. Therefore, the improvement of their thermal conductivity is a key factor for their exploitation in similar applications. Summer thermal performance of a DSHP coupled with a flat-panel GHE and a mixture of enhanced PCMs as backfill material is better than in wintertime, according to the higher temperature difference between undisturbed ground and GHE.
Therefore, DSHP coupled with a shallow GHE backfilled with PCMs having high thermal conductivity is an interesting alternative for the air conditioning of buildings, especially in high latitudes of Southern European countries, where cooling requirements are relevant. This solution may also avoid frosting problems and save the system against extreme temperatures that affect the air-mode, by using the ground as an alternative heat source, which makes the HVAC system more resilient against the incontrovertible climate change.
A future development of this research will be the experimental validation of the main conclusions and the assessment of the economic viability of the enhanced system. Author Contributions: All authors contributed equally to this work. All authors designed the simulations, discussed the results and implications and commented on the manuscript at all stages. M.B. conceptualized the study and methodology. F.J.G.G. carried out the formal analysis and led the development of the paper. Both authors performed numerical simulations and result analysis and discussion. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

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

Nomenclature c P
Specific heat capacity, J/(kg · K) c P eq Mixture specific heat capacity, J/(kg · K) c P i Specific heat capacity of component i, J/(kg · K) c P s Specific heat capacity of wet sand, J/(kg · K) c s