Economic Optimal HVAC Design for Hybrid GEOTABS Buildings and CO2 Emissions Analysis

Abstract: In the early design phase of a building, the task of the Heating, Ventilation and Air Conditioning (HVAC) engineer is to propose an appropriate HVAC system for a given building. This system should provide thermal comfort to the building occupants at all time, meet the building owner’s specific requirements, and have minimal investment, running, maintenance and replacement costs (i.e., the total cost) and energy use or environmental impact. Calculating these different aspects is highly time-consuming and the HVAC engineer will therefore only be able to compare a (very) limited number of alternatives leading to suboptimal designs. This study presents therefore a Python tool that automates the generation of all possible scenarios for given thermal power profiles and energy load and a given database of HVAC components. The tool sizes each scenario properly, computes its present total cost (PC) and the total CO2 emissions associated with the building energy use. Finally, the different scenarios can be searched and classified to pick the most appropriate scenario. The tool uses static calculations based on standards, manufacturer data and basic assumptions similar to those made by engineers in the early design phase. The current version of the tool is further focused on hybrid GEOTABS building, which combines a GEOthermal heat pump with a Thermally Activated System (TABS). It should further be noted that the tool optimizes the HVAC system but not the building envelope, while, ideally, both should be simultaneously optimized.


Introduction
In the early design phase, a design engineer will typically firstly estimate the energy profile necessary to heat and cool the future building based on a few parameters such as the building conditioned volume, its floor surface area, the type of occupancy, its windows-to-wall ratio, whether the structure is light or heavy, and typical weather conditions at the given location. Based on the energy profile and optionally some additional constraints set by the client, the design engineer will then propose a combination of properly sized heat and cold production and emission systems (further referred to as the Heating, Ventilation and Air-Conditioning (HVAC) scenario) able to deliver the required powers and energy loads. Finally, the investment, running, maintenance and replacement costs can be estimated along with the total CO 2 emission or another indicator for environmental impact. Ideally, several alternative HVAC scenarios are compared with each other to choose the most appropriate one. However, since the presented work flow is highly time-consuming, the number of considered alternative scenarios remains (very) limited leading to suboptimal designs. This work presents therefore a Python tool that automates the generation of all possible scenarios for given thermal power profiles and energy load and a given database of HVAC components. The tool further sizes each scenario properly, computes its present cost (PC) based on the investment, running, maintenance and replacement costs, and the total CO 2 emission associated to the building energy use. Finally, the different scenarios can be searched and classified to pick to the most appropriate scenario. The tool is based on static calculations based on standards, manufacturer data and basic assumptions similar to those made by engineers in the early design phase. It should further be noted that the tool optimizes the HVAC system but not the building envelope while, ideally, both should be simultaneously optimized.
This paper is structured as follows. Firstly, Section 2 describes common methods used to improve or optimize HVAC systems of buildings. Secondly, Section 3 introduces the methodology proposed in this work, Section 4 clarifies its assumptions, and Section 5 lists the considered HVAC components and summarizes their characteristics. Finally, Section 6 compares the Present Cost (PC) and the CO 2 emissions of conventional, (pure) GEOTABS, and economic optimal HVAC designs for different boundary conditions, and Section 7 concludes and proposes some further improvements.

Existing Methods
In Europe, buildings are responsible for 40% of the total energy use from which half is used for heating and cooling [1]. In accordance to the European Union's Directive 2010/31/EN [2], the energy requirements for buildings are becoming increasingly stringent. Not only should the buildings become more energetically efficient through a better design of the building envelope (insulation, shading, etc.) and the use of more efficient HVAC devices, but buildings are also obliged to use a minimum of renewable energy sources. Both requirements put building designers and installers under pressure to be innovative and up-to-date with the various and rapidly improving available technologies. Furthermore, systems are becoming increasingly complex. Where a standard gas-or oil-fired boiler and a compression cooling machine connected to radiators and ventilation units used to be installed everywhere, buildings using a combination of production systems such as gas boilers, heat pumps, etc., emission systems such as radiators, fan coil units, chilled beams, TABS, etc., and additional renewable energy sources such as photovoltaic panels, solar boilers, borefields, etc. are becoming common in Europe. With this increase in system complexity and the lack of knowledge about these new and rapidly improving technologies, design tools are becoming highly necessary to help designers and installers to optimize both the design and the control of buildings. According to Ellis and Mathews (2002) [3], researchers believe that design tools based on an integrated approach where the system efficiency is optimized taking the interactions between the various factors and their constraints into account could lead to savings around 70%.
According to the American Society of Heating, Refrigerating and Air-Conditioning Engineers (ASHRAE), three types of equipment related programs exist to help HVAC designers: (i) electronic catalogues that simply list available components respecting given constraints, (ii) equipment optimizers that propose a range of possible equipment alternatives for given design criteria, and (iii) equipment simulation programs that can take the interactions and constraints of different devices connected to each other into account [3]. While only equipment simulation programs can truly help the design engineer to optimize the system, such a tool is confronted to an intricate problem. Firstly, the number of design variables (type, size and combination of systems) is usually very high. Secondly, the operation efficiency of each system is very dependent on its components and how they interact, and, finally, the design is heavily based on the building's energy demand, which is difficult to predict in the early design phase.
In the literature, building and HVAC designs are typically optimized by heuristic methods such as genetic algorithms. Due to their slow convergence, generally only a limited number of parameters are optimized while the control is fixed or optimized over a short period of time. Another popular approach is to formulate the optimal design problem as a mixed integer linear programming (MILP) problem, allowing many variables to be solved in a reasonable computing time. However, a MILP formulation requires linear models that prevents its use in simulation based building design optimization. Nguyen et al. [4] provide an extensive literature study on building optimization. Examples of optimization using genetic algorithms can be found in e.g., Ooka and Komamura (2009) [5], who used a genetic algorithm to optimize the equipment capacity and operation planning of buildings, and, in Seo et al. (2014) [6], who used a multi-island genetic algorithm to optimize the HVAC system of apartments. The MILP approach has been used by, e.g., Ashouri et al. [7] to simultaneously optimize HVAC equipment, sizing and operation, and by Patteeuw and Helsen [8], who developed a similar method but also included multiple temperature levels to represent energy storage and conversion efficiencies in a more realistic way, and the explicit modelling of the electricity generation side. Heuristic and MILP optimization methods can also be combined. Evins [9] simultaneously optimized the building envelope of an office building, its HVAC, and its control, using a multi-objective, multi-level optimization scheme. At the top level, a genetic algorithm optimizes some building envelope parameters and the choice and size of the HVAC system by running an EnergyPlus model. At the lower level, a MILP algorithm is used to optimally control the HVAC with minimization of the running cost as objective. Due to the huge problem formulation, the author used a combination of a workstation with twin eight-core processors and 64 Gb RAM and additional 960 computing cores each with 4 Gb RAM to solve the MILP problems in parallel.
On the more practical and political sides, EU's Directive 2010/31/EN compels the EU countries to limit the energy used by the building sector but leaves each country the freedom to develop its own legislative framework to reach the energy goals. In Belgium, a software called EPB (in Dutch and PEB in French) (version 9.0.0, Vlaamse energieagentschap (VEA), Belgium) is used to compute an energy level indicator. The energy level is evaluated based on quasi static energy loss computations and on a point table to grade different HVAC systems. The Flemisch Energy Agency (VEA) publishes every second year a report that investigates for which energy level the building would also be cost optimal with respect to the sum of the investment, running and replacement costs [10]. While the results are very instructive for decision makers of Flanders, the method is not practical for designers as it is difficult to relate the energy indicator to real energy use and the designer is restricted to the initial assumptions made by the program. Furthermore, the EPB software is not able to calculate dynamic aspects, system integration and control correctly.
The drawback of the presented methods is that they require intensive computer power, specialized software, and intensive modelling and set up work. The results are further difficult to analyse as they are influenced by numerous optimization parameters and assumptions/approximations made by the tools. Therefore, this work proposes to go back to the straightforward, but time-consuming method used by engineering offices and to automate it. The method should be based on a limited set of parameters, and it should remain intuitive to be usable in the early design phase. The starting point of the method is the heating and the cooling load duration curves (LDC) of the building, which represent the thermal powers required to condition the building. The thermal powers are furthermore ordered and each power is plotted with a width equal to the sum of hours that it is used, such that the integral of the LDC represents the yearly energy load of the building. The developed method proposes then all possible HVAC scenarios composed of the devices present in the database, which can provide the powers and the energy load of the LDC and computes the scenarios PC (including the investment, running, maintenance and replacement costs) and CO 2 emissions (see Section 3). The method has the main advantage to automatically generate and size, and to return the necessary information about each possible HVAC scenario such that design engineers can make the optimal choice, while the method complexity is sufficiently low such that the results can be easily manually verified.

Methodology
The methodology used by the developed tool to optimize the HVAC system for a given building is schematically represented in Figure 1, which shows the four steps that are consecutively carried out.
Step 1 Step 2 Step 3 Step 4 [W] [h] Energy allocation Adapt energy Figure 1. Schematic representation of the optimal Heating, Ventilation and Air-Conditioning (HVAC) design methodology split in four consecutive steps. The generation of load duration curves (LDC) based on a limited set of building parameters is not yet implemented and is therefore represented as a dotted grey arrow.
Step 1: Gather building and components data Firstly, a database of HVAC components is created containing for each component (see Section 5): its cost function returning its price as a function of its size (kW, m 2 , etc.), its maintenance cost, and its life expectancy (see Table 1 and Figure 3). If the component is a production device, the energy efficiency is given for different temperatures (for heat production: 30, 45, 60 • C, and for cold production: 10, 15, 20 • C, see Table 2) and if the component is an emission device, its nominal supply and return temperatures are specified (see Table 3). The database also contains the component constraints (e.g., ground source passive cooling (GSPC) cannot deliver water colder than 15 • C). Additionally, a number of building parameters are required such as its total floor area (used to size the Concrete Core Activation (CCA) system) and the ground characteristics that influence the borefield design. In its final state, the tool should also be able to estimate the energy profile of the building based on its parameters in the form of load duration curves (LDC). However, this is not yet implemented and the building's LDCs are thus now assumed to be given.
Step 2: Generate and size all possible HVAC scenarios The second step is to generate all possible HVAC scenarios composed of the devices present in the database, which can deliver the energy and thermal powers of the given LDC using the HVAC components of the database. In order to limit the number of possible scenarios, a scenario can contain a maximum of two heat and two cold production systems, each coupled to its own emission system. Note that reversible components count for two: one for heating, one for cooling. Figure 2 illustrates how the different scenarios are generated. Firstly, a set of production and emission systems is taken from the available components of the database. Secondly, the production systems are sorted by decreasing efficiency and the heat (cold) emission system with the lowest (highest) supply temperature is coupled to the most efficient heat (cold) production system. The maximum heating (cooling) power present in the LDC can now be divided between the two heat (cold) systems according to predefined hybrid fractions (currently set as [0, 25, 50, 75, 100]%). Note that both the heating and the cooling hybrid fractions can vary independently, that coupled production and emission systems are sized together, and that, when a 0% hybrid fraction is allocated, this means that the component is not installed. When a reversible system is installed, it will be sized by the maximum of its heating and cooling powers, taking its cooling to heating power conversion factor (see Table 3) into account. This is necessary as, for example, the cooling power of a high temperature 4 pipes fan coil unit (FC4P-HT) at its nominal conditions is only 60% of its heating power due to the smaller temperature difference between the air of the zone and the fan coil supply water for cooling than for heating.  While most of the components are sized according to their nominal power, CCA and Floor Heating (FH) are sized to a realistic fixed 80% of the building floor surface area and the ground heat exchanger (i.e., the borefield) is automatically designed according to the standard ISSO73 [32]. Standard ISSO73 returns the number of boreholes and their depth necessary to deliver a given maximum thermal power and a given amount of yearly energy and it takes into account the regeneration level, the ground and grout characteristics, the borefield configuration, the borehole type, the maximum borehole depth, and the minimum allowed ground temperature (T min ). The standard was originally developed to size heating dominated borefields based on the results returned by the Energy Earth Design (EED) software [33] with an accuracy of ±10%. In this work, the work flow proposed by the standard has been translated in the Python tool to automatically size the borefield and the method is extended as follows to allow cooling dominated cases. When the borefield is heating dominated, its sizing depends on the temperature difference between the initial ground temperature (T 0 ) and T min . In Belgium, this temperature difference is typically around 10 • C as T 0 is typically around 10 • C and T min is here chosen equal to 0 • C to avoid freezing of the ground. When the borefield is cooling dominated, the borefield is sized such that the ground temperature does not exceed a maximum value (T max ), which is chosen in this work to be 17 • C, to ensure that passive cooling with the ground can always be used. The same work flow can then be applied to size cooling dominated borefields by simply using the cooling energy and power instead of the heating ones and by choosing T min so that T 0 − T min = T max − T 0 , which means T min = 3 • C for the mentioned assumptions.
Finally, the original energy of the LDC, which corresponds to the building energy needs for space heating and cooling, is corrected to include the emission system efficiency as specified by the standard NBN EN15316-2-1:2007 [34] (see Section 5 and Table 3).
Step 3: Compute energy use, costs and emissions Once the type and the size of each HVAC component of a scenario are defined, the energy from the LDC can be allocated to the different production components such that the most efficient component delivers the base load. Notice that Figure 1 shows only one LDC while actually two LDCs are needed: one for heating, one for cooling. The scenario present cost (PC) can now be computed as: with I 0 the initial investment cost calculated as the sum of the different component cost functions (see Table 1), N the number of years during which the system will be used (in this work, 20 years is assumed), C e , C m , and C r the energy, maintenance, and replacement costs computed using the efficiencies, maintenance percentages, and life expectancies of Tables 1-3 and a given electricity and gas price of 0.1466 and 0.0534 e/kWh [35,36], and V r the remaining value of the system after N years. Note that the investment cost of the components and the energy prices are assumed to stay constant over the years and no inflation is considered. This assumptions are chosen due to the too high uncertainty on those factors and in order to keep the computations more transparent. However, the different costs are discounted with discount factor R j = 1 (1+i) j with j the number of years counted from the investment time and i the discount rate. An interest rate of 4% is used in this work. V r is further calculated linearly with the number of years the device can still be used except for radiators, CCA and the ground heat exchanger (GHEX-V) for which V r is zero as they cannot be easily dismantled and re-used.
Besides the PC, the total CO 2 emissions caused by the energy use of the building are computed using the primary energy conversion factor of 0.056 kg/MJ for gas and 0.179 kg/MJ for electricity [10].
Step 4: Choose the most appropriate HVAC system The final step compares the PC and the total CO 2 emission of each scenario over a realistic 20 years lifetime of the building such that the design engineer can make the most appropriate choice.

Assumptions and Limitations
The main assumption/limitation of the presented methodology is that the LDC of the building is assumed given and therefore no dynamic simulations have been carried out. In the future, the method should include an LDC generator. The powers and energy of the LDC furthermore depend on the chosen HVAC scenario, which is taken into account by the emission efficiencies from NBN EN 15316-2-1:2007 [34]. However, the accuracy of such efficiencies should be investigated by dynamic simulations. Finally, as no dynamic simulations are carried out, comfort cannot be explicitly checked. The method assumes that, if the HVAC system can provide the thermal load defined by the LDC, comfort will be guaranteed and it is up to the designer to choose the most comfortable HVAC system, while taking into account the energy use, costs and emissions returned by the method.
The following assumptions are also made: • While different efficiencies for different supply temperatures are used, the nominal power of production devices does not depend on their supply temperature. This might cause an over or underestimation of the device investment cost when it is used in different nominal conditions. • For geothermal systems, only the most common vertical ground heat exchangers (GHEX-V) are considered. Furthermore, it is assumed that Ground Source Heat Pumps (GSHP) are not reversible and that cooling can only be done passively using a GSPC. This assumption might cause an overestimation of the investment cost of GSHP systems.

•
If CCA or FH is installed, it is assumed that it will be installed for the whole building and cover a fixed 80% of the building floor surface area. • So far, the ventilation system and the building climate controller are not included in the HVAC design optimization while they may have a significant impact on the investment cost, energy use and achieved comfort. In addition, the number of considered HVAC components is limited (see Table 1).

•
No inflation is taken into account due to its large uncertainty, but all costs are discounted as described in Section 3.

Components Information
As described in Section 3, the first step to optimize the HVAC system is to gather the necessary information about its components. The nomenclature, the cost function, the maintenance cost, and the life expectancy of each considered HVAC component are listed in Table 1 and the cost functions are also represented in Figure 3. The cost functions are derived from the (installed) cost data of 18 recently built buildings from two different Belgian engineering offices. The maintenance costs, expressed as a percentage of the investment cost, and the life expectancies are retrieved from the standard NBN EN 15459 [37]. The database can easily be extended with new components and the required parameters are known.   Table 1. Table 2 contains the assumed energy conversion efficiencies of the different heat and cold production components for different supply temperatures, based on the average of values found in technical sheets of different manufacturers and in the Annex 48 report [39]. When the efficiency is not given, it means that the component cannot supply water at that particular temperature. For the emission components, Table 3 indicates their nominal conditions for heating and cooling, the ratio between their cooling and their heating power and their efficiency. The emission system effiencies are taken from NBN EN15316-2-1:2007, and they represent the system additional heat losses. These heat losses are due to the non-uniform internal temperature distribution in the conditioned zones (caused by stratification, heat emitters along outside wall/window, differences between air temperature and mean radiant temperature), the non-ideality of the operative temperature control (causing temperature variations and drift) and the extra heat losses of heat emitters embedded in the building structure towards the outside. Table 2. Efficiencies (i.e the ratio between the useful thermal energy and the used electrical energy) of production components based on values from technical datasheets (TS) and Annex 48 [39]. In the case of Ground Source Heat Pump (GSHP), the given efficiency corresponds to the Seasonal Performance Factor (SPF) assuming a circulation pump power of 2.5% of the heat pump thermal power. The Coefficient of Performance (COP) is given in brackets and it is evaluated for a source temperature of 5 • C.   Table 1). Nominal conditions are given as T s /T r /T a with T s , T r and T a the supply, return and zone air temperatures. The power ratio column indicates the ratio between cooling and heating powers for reversible systems and the last column introduces a correction factor (η HVAC ) for the energy use of the emission system to represent the suboptimality of its control.

Results
In this section, two types of scenarios are considered: (i) the so-called minimum scenarios (M-Sce) that are composed of one single production system for heating and one for cooling, each coupled to its own emission system (see Section 6.1), and (ii) so-called hybrid scenarios (H-Sce) composed of two heating and two cooling systems (see Section 6.2). For both M-Sce and H-Sce, their PC over 20 years and their yearly CO 2 emissions are compared for different parameter values. In the case of M-Sce, the PC only depends on the maximum heating (P h,max ) and cooling (P c,max ) powers, and on the total yearly amount of heating (E h ) and cooling energy (E h ) delivered to the building. Additionally, the building floor surface area (A floor ) and the ground thermal conductivity (λ g ) can influence the scenarios PC if it uses CCA and a GSHP. In the case of H-Sce, the shape of the LDC is also important as it defines the ratio between the energy covered by the base load system and the energy covered by the peak load system.
In order to limit the degrees of freedom, the parametric space is represented in this work by five parameters: P h,max , FLH h , α, P spe , and λ g . FLH h is the number of full load hours that the heat production system runs per year. Multiplying FLH h with P h,max gives thus E h . The parameter α is defined as α = E c /E h and its value indicates whether the building is heating (α < 100%) or cooling dominated (α > 100%) or in thermal balance (α = 100%). Note that a GEOTABS-scenario (GEOTABS-M-Sce) (i.e., a scenario only composed of a GSHP, a GSPC and CCA) with α = 100% does not have a borefield with a balanced heating and cooling load as a fraction of the heat comes from the electricity use of the GSHP and not from the ground. In this work, α is varied by changing the ratio P c,max /P h,max while the number of full load hours of the cold production device (FLH c ) is always kept equal to FLH h . The influence of A floor is represented by the maximum specific power per unit floor area: P spe,h = P h,max /A floor . Finally, λ g is used as an indicator of the ground suitability for GSHP systems. The values of these five parameters are taken from realistic ranges for large Belgian buildings (see Table 4). The range for P h,max is further chosen such that the cost functions of the production components are always evaluated within their validity ranges, even when hybrid systems are installed (see Section 6.2), and P spe,h is kept equal to 27 W/m 2 for M-Sce to ensure that the maximum heating and cooling powers can be provided by CCA without auxiliary emission systems, even for the maximum value of α = 150%.
For the reader's convenience, the meaning of the different symbols and abbreviations used thoroughout this chapter are summarized in Tables 1, 4 and 5. Figures 5-10 furthermore use the following convention for the legend entries: each keyword separated by '/' represents an HVAC component installed in the scenario, and the number following the keyword indicates its nominal power as a fraction of P h,max for heating devices and P c,max for cooling devices. For example, for a LDC with P h,max = P c,max = 150 kW, GSHP0.75/CGB0.25/GSPC1/R-HT0.25/CCA means that the scenario is composed of a GSHP of 112.5 kW, a CGB of 37.5 kW, R-HT with a total power of 37.5 kW and CCA (for which no fraction is given as CCA is assumed to always cover the whole building floor surface area). (*1) this value is only used in Section 6.2, not in Section 6.1. Table 5. Description of the abbreviations used in Section 6. When the abbreviation is followed by a star (*), it indicates that the scenario is cost optimal compared to the scenarios of the same type.

Abbr. Description
Sce HVAC scenario(s) containing some heat/cold production and emission systems.

M-Sce
Minimal Sce: Sce containing only one heat, one cold production system, each coupled to its own emission system. (*1)

H-Sce
Hybrid Sce: Sce containing max two heat and max two cold production systems, each coupled to its own emission system. (*1) GEOTABS-M-Sce M-Sce only composed of GSHP/GSPC/CCA.

Geo-Sce
Sce with at least one geothermal production system.
(*1) In case of reversible emission system, one heat and one cold production system can be coupled to the same emission system.

Results for Single Production/Single Emission Scenarios
In this section, scenarios (M-Sce) composed of a single heat and single cold production system, each coupled to its own emission system, are considered. Figure 4 plots the PC of each M-Sce for the different parameter values of Table 4. Each row assumes a different FLH h and λ g , each column a different α while the x-axes specify P h,max , and each line represents a different M-Sce. The lines representing GEOTABS-M-Sce are plotted in green with square markers and the lines referring to the most conventional scenario type (Conv-M-Sce) (composed of a CGB, a CCM and FC4P-HT) are plotted in red with diamond markers. Finally, the background of the graphs is set in blue if GEOTABS-M-Sce extract more heat from the ground than they inject (i.e., the borefield is heating dominated) and in orange otherwise.
From Figure 4, it appears that, in general, Conv-M-Sce are the most cost efficient M-Sce that can be installed. When only a low cooling load is required (α = 50%), Conv-M-Sce are a factor 1.0 to 1.7 cheaper than GEOTABS-M-Sce and for cooling dominated buildings (α = 150%), the factor varies between 1.2 to 1.8. Only the cases when the cooling and heating loads of the borefield are in balance (i.e., α = 82% for a COP = 5.6) and when the building uses a minimum of 1750 FLHs, GEOTABS-M-Sce are financially advantageous (up to a factor 1.1 for the highest λ g and FLH). Interestingly, the lines representing the different M-Sce do not cross each other for different P spe,h , which means that the cost optimal M-Sce (M-Sce*) for the investigated parameter values do not depend on P spe,h . The cooling dominated columns of Figure 4 (α ≥ 100%) further indicate that other M-Sce than Conv-M-Sce can be cost optimal. This is illustrated in Figure 5.    Table 4. The blue backgrounds indicate that borefield of the Geo-Sce is heating dominated (more energy extracted than injected). The lines representing GEOTABS-M-Sce are plotted in green with square markers and the lines referring to the most conventional scenario type (Conv-M-Sce) are plotted in red with diamond markers. Figure 5 contains the same information as Figure 4 evaluated at P h,max = 150 kW, but the CO 2 emissions are added to the figure. The left y-axis represents the PC (filled markers) of the M-Sce (in e/20y) and the right y-axis returns the CO 2 emissions (unfilled markers) in ton/1y, the x-axes indicate the number of FLH h and λ g , and each subplot corresponds to a different α value. For each condition, the GEOTABS-M-Sce and the Conv-M-Sce are plotted with green square markers and red diamond markers, respectively. Additionally, M-Sce* and the cost optimal geothermal scenario (Geo-Sce*) are also included if they differ from GEOTABS-M-Sce and Conv-M-Sce. Finally, the CO 2 emissions of each plotted scenario are represented by the same (but unfilled) marker, and their value is to be read on the right axis.   Figure 5. Present Cost (PC , filled markers) and CO 2 emissions (unfilled markers) of HVAC scenarios (Geo-Sce) containing only one heat and one cold production system, each coupled to its own emission system for the different parameter values of Table 4 and P h,max = 150 kW. Only the GEOTABS (GEOTABS-M-Sce), the most conventional (Conv-M-Sce), the cheapest (M-Sce*) and the cheapest geothermal scenarios (Geo-Sce*) are plotted. The legend specifies the components of the HVAC scenario and the fraction of installed power they cover. Figure 5 shows that, when a large amount of cooling is required (α is high), GEOTABS-M-Sce are systematically (much) more expensive than Conv-M-Sce despite the fact that cooling with GEOTABS-M-Sce is more energy efficient (EER = 20) than with Conv-M-Sce (EER = 3.9). This is due to the fact that the borefields of GEOTABS-M-Sce with α ≥ 100% are sized according to the cooling load, which leads to larger, thus (much) more expensive, borefields. By sizing the borefield only for the heating load and covering the cooling load with a CCM (GSHP1/CCM1/CCA1 markers in Figure 5), the scenario price could be decreased by a factor up to 1.3. The cost of the Geo-Sce can further be reduced by allowing a hybrid design such that the borefield would deliver only a fraction of the cooling load (see Section 6.2). Similarly, the FC4P-HT are becoming very expensive when they are used for unbalanced loads. The scenarios GSHP1/CCM1/CCA1 and those combining GSHP, GSHPC, R-HT, and FC2P-coo-HT can therefore become cheaper as the use of the expensive FC4P-HT can then be avoided. In the case of a thermally balanced borefield (α = 80%), the GEOTABS-M-Sce become very advantageous and their PC drops below those of Conv-M-Sce for buildings using at least 1750 FLH while they emit 58% less CO 2 . Finally, the use of R-LT and FC4P-LT instead of their high temperature equivalent appears not to be advantageous to be coupled to the condensing gas boiler and cooling machine as their additional investment costs are not compensated by the realized energy savings.

Results for Hybrid Scenarios
This section is similar to Section 6.1, but hybrid scenarios (H-Sce) are now allowed with a maximum of two heat and two cold production systems, each coupled to its own emission system (see Section 3). H-Sce depend on the same parameters as M-Sce (see Table 4), but the shape of the LDC now also influences their PC and CO 2 emissions. Therefore, two types of LDC are used, i.e., triangular and logarithmic, as represented in Figure 6. Furthermore, two different values of P spe,h are considered in this section (27 and 70 W/m 2 , see Table 4), as the CCA can now be complemented with a second emission system both for heating and cooling.  The following figures representing the PCs and CO 2 emissions of H-Sce for the different boundary conditions contain the following information. Firstly, as H-Sce can become a M-Sce if only 0% and 100% hybrid factions are used, green squares and red diamonds still represent GEOTABS-M-Sce and Conv-M-Sce, respectively. Note, however, that, when P spe,h = 70 W/m 2 , H-Sce will not contain any GEOTABS-M-Sce as the CCA cannot provide the thermal power without an auxiliary emission system. Secondly, the cost optimal geothermal scenario (Geo-Sce*) is also plotted in each figure. Finally, if the overall cost optimal H-Sce (H-Sce*) does not coincide with an already plotted scenario, H-Sce* is added to the figure.

Triangular LDCs
The shapes of the triangular LDCs are chosen such that they correspond to the P h,max , FLH h , and α values of Table 4. Figures 7 and 8 assume a P spe of 27 W/m 2 and 70 W/m 2 , respectively, resulting in different optimal scenarios as the CCA (if installed) is a factor 2.6 cheaper for the latter (CCA are assumed to cover a fixed 80% of the building floor surface area), but it needs to be complemented with an auxiliary emission system to provide enough thermal power.  Figure 7. Present Cost (PC , filled markers) and CO 2 emissions (unfilled markers) of HVAC scenarios (Geo-Sce) containing maximum two heat and two cold production systems, each coupled to its own emission system for the different parameter values of Table 4 but P h,max = 150 kW, P spe = 27 W/m 2 and the Load Duration Curve (LDC) are of triangular shape. Only the GEOTABS (GEOTABS-M-Sce), the most conventional (Conv-M-Sce), the cheapest (M-Sce*) and the cheapest geothermal scenarios (Geo-Sce*) are plotted. The legend specifies the components of the HVAC scenario and the fraction of installed power they cover.
In the case when P spe = 27 W/m 2 , the results for Conv-M-Sce and GEOTABS-M-Sce are the same as in Section 6.1. Figure 7 shows that, for heating dominated buildings, Geo-Sce* are generally composed of a combination of a GSHP (25%) and a GSPC (25%) coupled to FC4P-LT and of a CGB (75%) and a CCM (75%) coupled to FC4P-HT. Here, the Geo-Sce* remain a factor 1.03 to 1.28 more expensive than the Conv-M-Sce but emit 20% less CO 2 than Conv-M-Sce. For cooling dominated cases, installing a GSPC covering only a fraction of the required cooling power (25 to 50%) significantly reduces the PC cost differences between Conv-M-Sce and Geo-Sce* to such an extent that Geo-Sce* become the H-Sce* for buildings with large energy loads (2500 FLH). Furthermore, Geo-Sce* systematically reduce the CO 2 emission by more than 50% compared to Conv-M-Sce.
As illustrated by Figure 8, in the case when P spe = 70 W/m 2 , Geo-Sce coupled to CCA and FC4P-HT are the optimal scenarios H-Sce* when α ≥ 80 and have similar cost as Conv-M-Sce for α = 50%. For example, Geo-Sce* can be up to 1.1 cheaper than Conv-M-Sce when α = 80% and up to 1.2 when α = 100% while Geo-Sce* save in both situations around 50% of the CO 2 emissions.

Logarithmic LDCs
LDCs have rather an (inverted) logarithmic shape than a triangular shape, and this paragraph therefore investigates the cost optimality of HVAC scenarios for LDCs described by f (x) = γ − β(x − ξ) with x the number of hours and f (x) the heating power. The parameter ξ influences the steepness of the curve and γ and β are calculated such that the different LDCs correspond to the parameter values of Table 4, and such that the LDCs tail is cut at P = 0.2P max (see Figure 6). The results are again analysed for P spe = 27 and 70 W/m 2 (see Figures 9 and 10).   Figure 8. Same as Figure 7 but P spe = 70 W/m 2 .  Figure 9. Present Cost (PC , filled markers) and CO 2 emissions (unfilled markers) of HVAC scenarios (Geo-Sce) containing maximum two heat and two cold production systems, each coupled to its own emission system for the different parameter values of Table 4 but P h,max = 150 kW, P spe = 27 W/m 2 and the LDC are of logarithmic shape. Only the GEOTABS (GEOTABS-M-Sce), the most conventional (Conv-M-Sce), the cheapest (M-Sce*) and the cheapest geothermal scenarios (Geo-Sce*) are plotted. The legend specifies the components of the HVAC scenario and the fraction of installed power they cover.
The results for logarithmic LDCs are very similar to those for the triangular LDCs. However, where the Geo-Sce* from the triangular LDCs would most of the time size its GSHP to 75% and its GSPC to 50 or 25%, the Geo-Sce* from the logarithmic LDCs would rather size them to 50% and 25%, respectively. This is a logical result as the number of FLH for the higher thermal powers decreases faster for logarithmic LDCs than for triangular ones. This also results in a lower PC for the Geo-Sce with logarithmic LDC while Conv-M-Sce and GEOTABS-M-Sce remain the same. For example, the Geo-Sce* is 1.25 times cheaper than Conv-M-Sce when FLH=2527 h, λ g = 2.4 W/(m K) and P spe = 70 W/m 2 and still saves 50% of the CO 2 emissions.  Figure 10. Same as Figure 9, but P spe = 70 W/m 2 .

Conclusions
This chapter presents a tool developed in Python to automatically compute the present cost (PC) and CO 2 emissions of all possible HVAC designs of a given building composed of the devices contained in the database. The tool was used to investigate when HVAC systems using a vertical ground heat exchanger (GHEX-V) become economically interesting and to compare their CO 2 emissions to those of a conventional HVAC system which uses a condensing gas boiler, a compression cooling machine and reversible high temperature four pipes fan coil units.
It was found that pure GEOTABS systems (i.e., an HVAC system composed of a GSHP system coupled to CCA) for buildings with a thermal power between 150 and 250 kW are generally 1.0 to 1.8 times more expensive than a conventional HVAC system when the GHEX-V is not in thermal balance but that pure GEOTABS systems have the lowest PC of all possible scenarios when ground thermal balance is reached and the building uses its GSHP system for at least 1750 full load hours (FLH) per year. For all cases, pure GEOTABS systems save more than 50% of CO 2 emissions compared to a conventional one.
It was further shown that, instead of installing pure GEOTABS systems, combining a GSHP system with TABS to cover the base load of the building and a condensing gas boiler and a cooling machine with fan coil units and/or radiators, results in hybrid GEOTABS systems with a PC that is lower than or similar to the PC of a conventional system, while still between 20% to more than 50% of CO 2 emissions are saved. Hybrid GEOTABS concepts can thus be very advantageous both economically and for the environment.