Comparison of different configurations of a solar-assisted ground-source CO 2 heat pump system for space and water heating using Taguchi-Grey Relational analysis

This work deals with the comparison of different optimized configurations of solar-assisted ground-source CO 2 heat pump systems using the Taguchi method and Grey Relational Analysis (GRA) when used for simultaneous space and water heating in cold coastal climate conditions. The configurations studied include: (1) the solar collectors (SCs) and borehole heat exchangers (BHEs) connected in series, with the working fluid flowing to the SCs first; (2) the SCs and BHEs connected in series, with the fluid flowing to the BHEs first; and (3) the SCs and BHEs connected in parallel to the heat pump. Eight parameters were considered for optimizing the systems ’ design, including the BHE length, BHE spacing, BHE number, SC area, BHE-SC mass flow rate, space heating return temperature, heat pump discharge pressure, and heat pump ’ s outlet temperature. The seasonal performance factor (SPF), levelized cost of heating (LCOH), and the estimated maximum annual ground temperature change (GTC) were chosen as performance indicators to evaluate system performance. The system model was developed using Modelica and 27 simulation runs for every configuration were implemented according to the L 27 (9 3 ) Taguchi orthogonal array. Single objective optimizations were first performed using the Taguchi method to determine the parameter combinations that would optimize the SPF, LCOH, and GTC, separately. After that, multi-objective optimization was performed using Taguchi-GRA to determine the control factor combination that would give the optimal overall performance when all output variables are considered simultaneously and given equal importance. When the performance indicators were considered separately, simulations show that configuration 2 gave the best SPF (4.025), configuration 3 gave the best LCOH (0.124 USD/kWh), and configuration 1 gave the best GTC (100.257%), respectively. Analysis of variance (ANOVA) showed that the SPF is most sensitive to the heat pump ’ s discharge pressure and outlet temperature; the LCOH to the BHE length, BHE number, and SC area; and the GTC to all the BHE and SC sizing variables. Multi-objective optimization showed that configuration 1 performs the best, giving a grey relational grade of 0.6875, which is equivalent to an SPF of 3.267, and LCOH of 0.155 USD/kWh, and GTC of 100.021%. BHE length was found to be the most influential parameter to overall performance, irrespective of the configuration.


Introduction
Heat pumps are identified as a central pillar in the effort to decarbonize space heating, domestic hot water (DHW) production, and cooling in buildings and industries [1,2].The European Union recognizes them as necessary technology to replace existing gas boilers and reduce the reliance on Russian natural gas [3].By absorbing lowtemperature thermal energy from their surroundings and raising it to levels suitable for utilization, they offer an alternative to traditional building energy supply, reducing reliance on high-grade forms of energy like electricity and fuels.
Most heat pumps currently use F-gases, specifically Hydrofluorocarbons (HFCs), as working fluid [4].HFCs were introduced as replacements for the ozone-killing Hydrochlorofluorocarbons (HCFCs) and Chlorofluorocarbons (CFCs) due to their comparable performance, efficiency, low toxicity, and non-flammability, all without damaging the ozone layer.However, certain HFCs possess a significantly high global warming potential (GWP), with some being approximately a thousand times more potent than CO 2 [5].Recent initiatives, like the EU's F-gases regulation (EC517/2014) and the Kigali amendment to the Montreal Protocol prescribe the phase-down of HFCs in the coming years [6,7].
Hydrofluoroolefins (HFOs) are being promoted by the chemical industry as a low-GWP HFC substitute.However, their reduced atmospheric lifetimes, while reducing GWP, can impact thermal and chemical stability, a strength of HFCs.HFOs are known to break down into trifluoroacetic acid (TFA), which enters the Earth through rainfall.TFA has no known environmental degradation pathways in water and accumulates and persists in environmental aqueous phases [8].Research is ongoing to fully understand the effects of TFA on the environment [9] but some studies show that it could be harmful to human health and aquatic life [10].As HFOs become more widely adopted, concerns about the accumulation of TFA become more prominent.
Fortunately, HFOs are not the only low-GWP alternative to HFCs.Natural working fluids, such as carbon dioxide (CO 2 ), ammonia, butane, isobutane, and propylene [11], have recently been gaining increasing attention as environmentally-safe working fluids.Among them, CO 2 (R744) appears to be a very promising option due to its zero ozone depletion potential, low global warming potential, non-toxicity, nonflammability, good thermodynamic properties, and affordability [12][13][14].The modern use of CO 2 in a trans-critical cycle was first proposed by Lorentzen [15] in 1990.Ongoing developments concerning trans-critical CO 2 heat pumps and their applications were comprehensively reviewed by Song et al [16].
Air source heat pumps (ASHPs) dominate the heat pump market, constituting about 60 % of global sales in 2021 [3].However, they typically exhibit poor low-temperature heating performance and frosting of the heat exchangers [17], given the temporal mismatch between the demand for heat and the temperature of the air, which serves as their heat source.Ground source heat pumps (GSHPs) offer a more efficient alternative for indoor climate control applications since they use the heat from the ground, which remains at a nearly constant temperature despite weather variations [18,19].
Numerous studies have already been conducted on CO 2 GSHP systems.Jin et al. [20] introduced a hybrid CO 2 GSHP system for space conditioning and DHW production in warm climates.It uses the ground and ambient air as heat sinks for cooling and relies solely on the ground for heating.They showed that this system could reach coefficients of performance (COPs) ranging from 3 to 5.5 with DHW at 65 • C. In a comparative study [21], they showed that an R410 ASHP exhibits better cooling COPs, except in temperatures above 30 • C, where the hybrid system's ground-based cooling excelled.However, for heating, the R410 ASHP consistently gave higher COPs.Additionally, they investigated the optimal control strategy for gas cooler pressure to reduce thermal imbalances in the ground e. Wu et al. [22] developed a 7-kW prototype and a model of a liquid-to-air CO 2 GSHP to investigate its performance in residential applications for both subcritical and trans-critical operations.Their simulations show that relative to those of R410 GSHPs', the cooling COPs of CO 2 GSHPs were higher for entering liquid temperatures (ELTs) below 20 • C, when the system was operated in the subcritical region, and were lower for ELTs above 20 • C when the system was operated in a trans-critical cycle.Wang et al. [23] showed that the CO 2 GSHP could reach COPs similar to that of R134a if applied to air conditioning systems and tap water heating, where high outlet water temperature is desired.Emmi et al. [24] studied a two-stage GSHP system that uses CO 2 (with an ejector) and R1234zeE as working fluids and found that the COPs and the Energy Efficiency Ratios (EERs) it exhibited do not differ from those of standard high-temperature doublestage heat pumps available in the market.Bellos and Tzivanidis [25] investigated the effects of different heat pump and BHE parameters on the COP of a CO 2 GSHP system.They showed that the most influential parameters are the temperature of the heater outlet, the number of BHEs, and the ground temperature.Other studies have focused on direct-expansion CO 2 GSHPs [26][27][28][29].
One of the primary challenges associated with the use of GSHPs is its high installation cost [18].Moreover, efficiently operating a single heat source heat pump system continuously can be difficult [17].Adding another heat source, such as solar thermal collectors (SCs), could help overcome these problems.SCs enable the use of shorter Borehole Heat Exchangers (BHEs) and provide more operational options.While there have been a few works on CO 2 solar-assisted ground source heat pumps (SAGSHPs), the number of studies and their scope is still very limited.Kim et al. [30] investigated how various operating parameters could affect the performance of a residential CO 2 SAGSHP for space heating.Choi et al. [31] performed simulations to compare an R22 and a CO 2 SAGSHP.They discovered that the R22 SAGSHP had a more stable performance and a higher heating capacity, but the contribution of SCs was more notable in the CO 2 hybrid heat pump system.Under basic operating conditions, the heating COPs of the R22 and CO 2 heat pumps were calculated to be 3.15 and 2.24, respectively.Both studies only focused on relatively low temperatures, given that the systems were only designed for space heating, and only considered a setup where the GSHP and the SCs were connected in parallel to a tank thermal energy storage (TTES).Sazon and Nikpey [32] conducted a sensitivity analysis, parametric study, and long-term performance simulation on a CO 2 SAGSHP system for simultaneous space and DHW heating.They utilized the system's seasonal performance factor (SPF), levelized cost of heating (LCOH), and the ground temperature change (GTC) it induces as performance indicators.The study emphasized the significant influence of the heat pump's discharge pressure and outlet temperature on the performance indicators.Additionally, they illustrated the importance of managing the GTC to ensure good long-term performance.However, they did not identify the combination of design parameters nor the system configuration that would optimize performance.
There are studies that have compared different configurations/operations of hybrid solar-geothermal heat pump systems, but they all seem to have focused on conventional subcritical heat pumps.Some studies have favored the BHEs and SCs to be connected in series.Si et al. [33] showed that a SAGSHP system where the working fluid flows through the SCs first before going to the BHEs outperforms a system that uses the SCs to store heat in a TTES in the daytime that will be used to restore soil temperature at night.When the fluid flows through the SCs first, surplus solar energy can be transferred to the ground to mitigate the degree of cold accumulation [34].In contrast, Razavi et al. [35] concluded that the system in which the fluid flows first to the BHEs before going to the SCs achieves the highest annual average COP.In this configuration, the fluid can be reheated by the SCs, which can help improve the COPs [34].Dai et al. [36] investigated six cases and observed that the SCs and BHEs connected in series performed better than when connected in parallel in terms of soil temperature recovery and COP.Nonetheless, other studies favored the BHEs and SCs to be connected in parallel.Lee et al. [37] noted that for heating-dominated applications, connecting the BHEs and SCs in parallel decreased the energy consumption of their system by 13.8 % compared to the serially connected SAGSHP.In another study [38], they showed that the parallel configuration demonstrated a 16.3 % greater heating capacity and 12.6 % higher COP at the optimal flow rate and ELT and ground temperatures of 19 • C and 14 • C, respectively.Configurations and operation modes of SAGSHPs affect the system performance considerably so the optimal match between solar and geothermal energy should be determined to run the system efficiently and sustainably [34].An effective way to compare different system configurations is to determine how they will perform individually when their designs are optimized for the same boundary conditions.However, SAGSHPs are quite complex systems, and optimizing them concerning multiple design parameters and performance indicators would entail heavy simulation workload.
The Taguchi method is useful for reducing the simulation runs needed to optimize systems.It employs an orthogonal array experimental design with a single analysis of variance (ANOVA) and utilizes the signal-to-noise (S/N) ratio to assess the parameter settings that minimize the sensitivity of system performance to the sources of variations [39,40].Previous studies have utilized the Taguchi method in optimizing GSHP systems [18,[41][42][43][44][45].Verma and Murugesan [46] applied the Taguchi method to a conventional SAGSHP system to optimize BHE length and the solar collector area.
The Taguchi method only focuses on the optimization of one performance indicator at a time.To facilitate multi-criteria optimization, where multiple performance indicators are considered simultaneously, some studies have combined the Taguchi method with Grey relational analysis (GRA).Several studies have utilized Grey relational-based Taguchi analysis for various energy systems, such as solar energy systems [47,48], absorption refrigeration [49], organic Rankine cycle [50,51], a micro-gas turbine system with seasonal storage [45], and industrial gas turbines [52].
The growing use of heat pumps should go hand in hand with the use of environmentally-friendly working fluids to mitigate potential issues associated with the use of synthetic refrigerants.However, to encourage wider use of this technology, there is a need for more information on their performance and operational characteristics.This paper aims to contribute to addressing this need by focusing on the design of a SAGSHP system that uses CO 2 as working fluid.The study involves the comparison of optimized CO 2 SAGSHP systems, considering different design parameters, system configurations, and multiple performance indicators, namely the SPF, LCOH, and the GTC that they induce.To reduce the computational load of multi-objective optimization, the Taguchi-Grey method was employed.This research contributes valuable insights to facilitate the broader adoption of environmentallyresponsible heat pump technology.

System description
The thermal energy system model utilized in this study was constructed using Modelica via the Dymola v2021x [53] graphical user interface (GUI).The Modelica Standard Library (MSL) v4.0.0 was used.It includes a CO 2 heat pump, SCs, BHEs, a TTES, and space and water heating demand information.The CO 2 HP model was developed using the Thermal Systems library v1.6.1 [54] and subsequently calibrated using experimental data.The BHEs were modeled using the MoBTES library v2.0 [55][56][57] modified to work with MSL v4.0.0.The SCs and TTES were both modeled with the Buildings library v9.0.0 [58].Simulations that span a year were run to capture seasonal variations in some parameters and boundary conditions.Three system configurations were investigated in this study: (1) Series A: the fluid flows through the SCs first before going to the BHEs, (2) Series B: the fluid flows through the BHEs first before going to the SCs, and (3) Parallel: the BHE and SC are connected in parallel to the heat pump and the fluid flow is equally divided among them.

The thermal energy system model
An earlier study gives details on the development, calibration, and validation of the different component models in the energy system, as well as the processing of the demand data [32].The thermal energy system model includes a CO 2 heat pump that implements a trans-critical cycle, SCs, BHEs, a TTES, and representative components for modeling space heating and DHW demands.
The SCs and BHEs serve as sources of low-temperature heat.Vertically discretized single U-tube BHEs were modeled using the thermal resistance and capacity modelling approach (TRCM) [59].The local heat transport model that links the BHE wall to the global volume implements a finite difference modelling (FDM) approach that divides the local volume into concentric ring elements.The heat flux to or from a volume element is defined on the global level and distributed among the elements of the local model.More details on the modelling of the BHEs are given in [55].The solar collector was modeled according to the EN12975 test standard [60].
The base value of the mass flow rate of the fluid (propylene glycol-water mixture) circulating through the SCs, the BHEs, and the water side of the heat pump's evaporator (0.4 kg/s) was estimated using the American Society of Heating, Refrigerating, and Air-Conditioning Engineers (ASHRAE) methodology for sizing BHEs [61].
The system incorporates a controller that regulates the rotational speed of the heat pump's compressor to achieve the desired outlet temperature of the hot fluid produced by the CO 2 heat pump.Additionally, an electric backup heater was positioned downstream of the heat pump to ensure that the fluid reached the target outlet temperature prior to going to the TTES.The TTES model uses the stratified storage tank model in the Buildings library [58] that implements several volumes that exchange heat among themselves and with the ambient via conductions.
Another controller manages the flow exiting the bottom of the TTES, guaranteeing that the temperature there remains above 50 • C.This is crucial to maintain the DHW distribution system's temperature above the proliferation range of Legionella bacteria (20---45 • C) [62].
Although not explicitly modeled, the heat required for space and DHW heating was supplied through a heat exchange process with the distribution system.Hot water was withdrawn from both the top and middle sections of the TTES to provide the energy needed for DHW production and space heating, respectively.The withdrawal rates were managed by controllers using the anticipated temperature of the fluid after undergoing heat exchange with the distribution system as set points.This represents a simplification of how the heat demand is met by the system.The cold return fluid, used for both space and DHW heating, as well as from the bottom of the TTES, is recirculated, mixed, and then reheated in the heat pump.The temperature of the fluid after DHW production was set at 10 • C, assuming a 3 • C pinch temperature relative to the assumed city water temperature of 7 • C. In this study, the system's comfort level was taken into account by ensuring that the return temperature of the fluid, after providing heat for space heating, remains at 25 • C or higher.The specific return temperature for space heating is influenced by the design and efficiency of the distribution system.It is expected to be at least as high as the temperature required for thermal comfort, which typically ranges from 15 to 21 • C during winter [63].Here, the space heating return temperature was given a base value of 30 • C. A lower return temperature implies a more efficient heat transfer to the distribution system; a higher return temperature implies a less efficient one.
The size of the CO 2 heat pump in this study was fixed at 6.5 kW since it was calibrated against a prototype unit having this capacity [64,65].Hourly thermal energy demand data for space and DHW heating, obtained from a school in Stavanger, Norway, were used as the reference of the demand inputs to the model.Since the demand from this facility exceeded the 6.5 kW heat pump capacity, it was adjusted to align with the heat pump's limitations.The raw data were first normalized by dividing all values by the measured maximum demand and then scaled to 3 kW and 3.5 kW for space and water heating, respectively.Peak demands, which comprise less than 1 % of the total annual demand, were filtered out for simplicity.As a result, the actual magnitudes of energy demand from the school were not used, but the demand patterns were retained.Space heating demand exhibits seasonal variations, peaking during the winter and during school hours, while decreasing in the summer.In contrast, water heating demand remains relatively consistent throughout the year.Fig. 1 shows the demand profile used in this study and provides typical weekday demand profiles for different seasons.Spring and autumn demand profiles closely resemble each other.Summer demand is notably lower, whereas winter demand surpasses all others in magnitude.Incident solar radiation data for Bergen, Norway [66], a city adjacent to Stavanger, were utilized due to the unavailability of precise weather data for Stavanger in existing databases.
The model uses the DASSL solver in Dymola to solve differential algebraic equations, which gives outputs at variable time steps to balance accuracy and efficiency.The tolerance was set to 0.0001.

The thermal energy system configurations
Shown in Fig. 2 are the schematics of the Modelica models of the three system configurations studied in the Dymola GUI.In the first configuration (Fig. 2a), the SCs and the BHEs are connected in series, with the fluid flowing to the SCs first before going to the BHEs (Series A).This allows the storage of excess solar energy into the ground.When solar irradiation is available, the SCs heat up the cold water-side fluid coming from the evaporator of the CO 2 heat pump.The solar-heated fluid is then directed to the BHEs, where it either extracts or injects energy, depending on its temperature relative to the ground.Through the BHEs, the ground could function as both a heat source and heat storage.Some studies on SAGSHP with a conventional sub-critical heat pump [33,36] favor this configuration.The second setup (Fig. 2b) also connects the SCs and BHEs in series but with the fluid flowing through the BHEs first before going to the SCs (Series B).In this case, the fluid can be reheated by the SCs and the heat pump can provide direct heating for buildings without starting up [34].Razavi et al. [35] noted this to be the most effective setup when used for space and DHW heating for a house located in Iran.The third configuration (Fig. 2c) connects the BHEs and the SCs in parallel to the heat pump.Some studies favor this setup for heating-dominated applications [37,38].When solar irradiation is available, half of the flow from the evaporator of the CO 2 heat pump gets heated up in the SCs while the other half in the BHEs.When it is unavailable, all the working fluid gets heated up with the BHEs.Depending on the availability and intensity of solar irradiation, this could be one way to reduce the load of the system to the ground.

Methods
The performance indicators considered in this work are the SPF, the LCOH, and the GTC induced by the system after the first year of utilizing it.Eight (8) parameters were considered as inputs during the optimization of the system, namely the BHE length, BHE spacing, BHE number, SC area, SC-BHE mass flow rate, the space heating return temperature, the CO 2 heat pump's discharge pressure, and the CO 2 heat pump's outlet temperature.The Taguchi method was implemented to first perform single-objective optimizations for each performance indicator for each system configuration.The contribution of each parameter to the value of every performance indicator was also determined.Multi-objective optimization was then performed by combining the Taguchi method with GRA.The overall performances of the optimized configurations were compared with one another.

Taguchi method
The Taguchi method is an experimental optimization technique that applies the standard orthogonal array, which gives the optimal number and set-up of the necessary trial runs.This allows the determination of the best level of each parameter to optimize a given response factor.The different phases of the Taguchi method include (1) problem formulation, (2) design of experiments, and (3) analysis of results.
Problem formulation entails the determination of the output variables, the objective functions, and the control parameters and their levels.In this work, the output variables or performance indicators include the SPF, LCOH, and the 1st year GTC.The objective is to keep the   SPF and GTC high while maintaining the LCOH low.
The orthogonal array was selected based on the number of parameters and levels.It specifies the number of trial runs needed to get sufficient information about the system.The minimum number of trial runs to be conducted can be determined by: where N Taguchi is the minimum number of trial runs, NV is the number of parameters, and J is the number of levels.In this study, there are eight (8) parameters studied, each having three (3) levels, which means that the minimum trial runs needed is 17.The smallest standard orthogonal array for three (3) factor levels that meets this minimum is the L 27 orthogonal array, which consists of 27 trial runs.A full factorial combination will require 6561 experiments (3 8 ) per system configuration, which would be too time consuming and complex.Applying the orthogonal array will provide a good estimate of information as with the full factorial combination with just 27 runs per configuration (total of 81).A higher resolution than this can be chosen to allow the detection of higher-order interactions.However, this would require more runs.The parameter combinations needed for the runs were determined using the design of experiments function in Minitab [67].
Analysis of the results entails the calculation of the signal-to-noise (S/N) ratio for each run and the implementation of ANOVA.The S/N ratio is a measure of robustness, which is used to identify parameters that reduce the process or product variability by minimizing the effects of uncontrollable factors.The S/N ratios of the SPF and the GTC were calculated using the "higher the better concept" (Equation ( 2)), while the S/N ratio of the LCOH was calculated using the "the lower the better" concept (Equation ( 3)).
where y i refers to the input parameter.ANOVA was carried out to assess the significance of each parameter.It estimates the relative importance of each parameter by calculating their respective percentage contribution to the performance indicator.The degree of freedom (DOF), sum of squares (SS), mean squares (MS), F-ratios, P-values, and percentage of contribution were included in the ANOVA tables.In this study, Minitab [67] was used to conduct ANOVA on the S/N ratios.

The performance indicators
The performance indicators chosen for this study are the SPF, the LCOH, and the GTC induced after the first year of utilizing the system.Since the study requires running several cases, it was decided to limit the simulation time to 1 year.Of course, more representative values of the performance indicators could have been obtained if the simulations covered the whole lifetime of the system.Nonetheless, as shown in [32], the SPF and LCOH of this system would not vary so much yearly if the overall change in ground temperature is minimized.This was assuming that the boundary conditions remained relatively consistent every year.
To calculate the SPF, the total energy delivered by the system to the demand for one (1) year was divided by the total energy utilized to run the system (Equation ( 4)).This represents the efficiency of the system.The system spends energy to run the compressors, the circulation pumps, and the back-up heater.A higher SPF indicates better efficiency.
where E SH is the energy delivered for space heating, E DHW is the energy delivered by the system for DHW production, and E cons is the total electricity consumption of the compressor, the pumps, and the back-up electric heater.
To calculate the LCOH, the operation cost and thermal energy delivered from a year of simulation were first assumed to be the yearly operation cost and energy generation of the system.This assumed annuity was then discounted back to its present value using the discount rate (r) and added to the total capital cost of the system.The sum was then divided by the discounted value of energy generation.
where A is the annual cost of operations, r is the discount rate, n is the lifetime of the system, C is the capital cost, and E is the yearly energy generation of the system.The LCOH signifies the cost of producing 1 kWh of heating, and it is affected by both the cost of the equipment and the efficiency of operations.This is a useful indicator of the cost efficiency of the system.It could also be used to see how the system compares with other available technology.A lower LCOH is desired.The data and assumptions used for cost calculations are summarized in Table 1.
The GTC was obtained from the BHE model.The model calculates the average ground temperature as it is utilized as a heat source or sink.Since the Taguchi method requires the output variables to be positive, the GTC in this study was defined to be the quotient of dividing the ground temperature after one year of simulation by the initial ground temperature.
where T final is the temperature of the ground after one year and T init is the initial average ground temperature.
A GTC value greater than 100 % means the ground temperature went up, while a value less than 100 % means that the ground temperature went down.In this study, changes to the ground temperature were kept minimal, but a higher-the-better concept was applied since the system was only intended for heating.When consistent boundary conditions are applied to a SAGSHP system, the yearly change in ground temperature is expected to reduce after the first year as the driving force for heat exchange between the circulating fluid and the ground gets smaller.The 1st-year GTCs can be assumed to represent the maximum expected temperature decline for a year of operation.This value was used as an indicator of how the system will perform long-term without simulating its whole lifetime.A GTC close to or greater than 100 % is desirable.

The parameters and their levels
The study considered eight (8) parameters, as outlined in Table 2. BHE length, BHE number, and SC area were selected due to their direct impact on the capital cost of the system.BHE spacing was included to address the thermal imbalance induced by the system to the ground [43].Note that the SC area also contributes to this aspect.SC-BHE mass flow, CO 2 heat pump discharge pressure, and the heat pump outlet temperatures were chosen since they are parameters that can be controlled or specified in a real-world system.Lastly, the space heating return temperature was included to indirectly represent the effects of the performance of the heat distribution system.
Levels are discrete values or conditions chosen for each parameter in a Taguchi experiment that are tested to determine their influence on an output variable.The number levels and their values are typically determined from experiments or simulations, considering the practical range of values that a parameter can take.Typically, parameters are given two (2) to three (3) levels.The number of levels affects the resolution of the experiment.A higher number of levels provides finer resolution, allowing for more precise identification of parameter effects and their interactions.However, specifying more levels will require more simulation runs.
In this work, the number and values of the levels for each parameter were specified by considering the results of the parametric and sensitivity studies done in an earlier study [32].Three levels were considered for each parameter (Table 2).For most parameters, level 2 adopted the base values; level 1 represented the lower bounds; and level 3 indicated the upper bounds.
In [32], BHE length was varied from 50 m to 200 m to see how it changes the values of the performance indicators.It was seen that the ground temperature can decrease substantially with utilization when the BHE length is kept short while the LCOH can significantly go up with longer BHEs.The BHE length of 100 m was chosen as the base value considering that it offered a good compromise between the GTC and LCOH values while keeping the SPF reasonable.The lower and upper bounds were then arbitrarily set to be − 50 m (50 m) and + 50 m (150 m) of the base value, respectively.
A similar approach was applied in choosing the values for the SC area levels.As shown in [32], the SPF remains relatively constant for different values of the SC area, but the GTC and the LCOH can both significantly change.The base value of six ( 6) collectors (13.92 m 2 ) was chosen as a compromise between the desired GTC and LCOH values.The lower and upper bounds were then set to be − 3 (3 collectors) and + 3 (9 collectors) of the base value, respectively.BHE spacing mainly affects the GTC.However, values of this parameter can be constrained by the availability of the installation area.
Sometimes a big area can be used but sometimes it can be quite confined to a limited space.In this study, it was assumed that the area that can be allotted for the BHEs is not very big, with a BHE spacing base value of m and an upper and a lower bound of − 2 m (3 m) and + 2 m (7 m) of the base value, respectively.
The study took measures to ensure that the value ranges for both the SC area and the BHE spacing were selected in a way that would keep the GTC close to one (1).This careful selection was aimed at maintaining a consistent long-term performance of the system.
In [32], it was demonstrated that optimal ranges exist for the SC-BHE mass flow rate and the discharge pressure of the CO 2 heat pump.For the SC-BHE mass flow, a base value of 0.4 kg/s was chosen since it lies within the range where the SPF is maximized, the LCOH is minimized, and the GTC is kept relatively high.A lower and an upper bound of half (0.2 kg/s) and double (0.8 kg/s) of this base value were then specified for levels 1 and 3, respectively.The lower bound was ensured to be higher than 0.1 kg/s since below this, the SPF starts to significantly drop down.
The optimal value of the discharge pressure depends on the operating conditions of the CO 2 heat pump.The design condition of the prototype unit used to calibrate the CO 2 heat pump model in this work is between 8.5 and 10 MPa.Since it was seen that the optimal discharge pressure that maximizes the SPF and minimizes the LCOH can occur within this, a base value of 9.0 MPa was chosen for this parameter.The lower and upper bounds of the CO 2 heat pump's design condition were then used as values of level 1 and level 3 of this parameter, respectively.
The outlet temperature of the heat pump was kept higher than 60 • C to ensure that the temperature in the DHW distribution system is kept higher than the proliferation temperature of Legionella (20---45 • C) [62].The values of level 2 and level 3 were then set to be + 5 (65 • C) and + 10 (70 • C) of this, respectively.This is in consideration of DHW distribution systems that may require higher temperatures.
The space return temperature from the facility where demand data were obtained was at around 30 • C and hence was used as the base value of this parameter.The lower and upper bounds were then set to be − (25 • C) and + 5 (35 • C) of the base value, respectively, to represent different efficiencies of the distribution system.

Grey relational analysis
GRA, together with the Taguchi method, was employed to perform multi-objective optimization.To perform this, the performance indicators for every experimental run were first normalized within the range of 0 to 1. Normalization was implemented in two ways: "the higher the better" and "the lower the better".The SPF and the GTC were normalized using the former while the LCOH values were normalized using the latter.For the "higher the better", the original array was normalized as follows [48]: For the "lower the better" concept, the original array was normalized by [48]: where y i (k) is the normalized value, max(x i 0 (k)) is the maximum value of the x i 0 (k) value and min(x i 0 (k)) is the minimum value of the x i 0 (k).After normalization, the Grey relational coefficients (GRC) were determined by the following equations [48]: Δ min = min ∀j∈i min ∀k ‖y 0 (k) − y i (k)‖ where Δ 0i (k)(i = 0, 1, 2, ⋯27; k = 1, 2, 3) is the deviation value between y 0 (k) and y i (k), φ is the distinguishing coefficient with φ ∈ [0,1], Δ max is the minimum value of Δ 0i , andΔ min is the minimum value of Δ 0i .The distinguishing coefficient value of 0.5 is typically utilized [47][48][49][50][51][52].Subsequently, the overall Grey relational grades (GRGs) were computed by: where w is the weight of the response.In this study, the contributions of all performance indicators to the overall GRG were assumed to be equal, i.e., w is equal to 1/3.The calculated GRGs were then used to perform a Taguchi process (S/N ratio calculation using "the higher the better" concept, ANOVA) to determine the parameter set-up that would give the optimal performance when multiple performance indicators are considered simultaneously.

Results and discussions
The values of the parameter levels considered here were based on earlier studies implemented on Configuration 1 (Series A).In [74], a BHE length of 82 m was calculated following the ASHRAE procedures for sizing BHEs.The ground was assumed to have the characteristics of Slate (density = 2760 kg/m 3 , specific heat = 920 J/kg-K, thermal conductivity = 2.1 W/m-K), one of the common rock types in some parts of Norway [75] while the thermal gradient was assumed to be 0.0125 K/m, similar to that of some wells drilled in Bergen, Norway [76].In [32], the BHE length was varied from 50 m to 200 m, while keeping the SC area at 13.92 m 2 (6 collectors).It was seen that a 50 m BHE can support the ~ 18500 kWh annual heating demand but ground temperature can decrease by 0.5 • C from the first year alone.This temperature decline will decrease with time if boundary conditions are kept unchanged, but it will affect the long-term capability of the system to provide thermal energy.Oversizing the BHE to 200 m can give slightly better SPFs and lower ground temperature decline, but it can entail LCOH values greater than the upper bound of the LCOH of typical solar thermal combi heat pump systems (0.206 USD/kWh) [84].
In [32], the effects of changing the values of other parameters, aside from BHE length, are given in more detail.This was considered in choosing the values of the parameter levels in this study, ensuring that a good balance among SPF, LCOH, and GTC can be obtained.Among all the parameters considered, the heat pump's discharge pressure and outlet temperature exhibited notable effects on all performance indicators.The presence of optimal SC-BHE circulation rate and discharge pressure were also noted.
The main objective of the study is to compare three different configurations of a CO 2 SAGSHP system when they are designed optimally in terms of three different performance indicators (SPF, LCOH, and GTC) for similar boundary conditions.The parameters considered include the BHE length (A), BHE spacing (B), BHE number (C), SC area (D), the SC-BHE mass flow (E), the space heating return temperature (F), the CO 2 heat pump discharge pressure (G), and the CO 2 heat pump outlet temperature (H) (Table 2).The Taguchi method was performed first to implement single objective optimizations concerning individual performance indicators (sections 4.1 and 4.2).Here, the effects of the different parameters on the different performance indicators for the three considered configurations were quantified.Afterward, the Taguchi-GRA was implemented for multi-objective optimization, with the performance indicators given equal importance (sections 4.3 and 4.4).In this part, the contributions of the parameters on the GRG, which represents the overall system performance, were calculated.

Single-objective optimizations: Taguchi method
As given in Table 2, eight (8) parameters, each having three (3) levels, were considered here.Following the Taguchi design concept, an L 27 orthogonal array was chosen.The structure of the array is given in Table 3.Each computational run was performed according to the combination of parameters determined by this array.The performance indicators obtained from simulating each of these runs for all three configurations are summarized in Table 4 while their corresponding S/N ratios are given in Table 5.The S/N ratios for the SPF and GTC were calculated using the "higher the better" concept while those for the LCOH were calculated with the "lower the better" concept.The S/N ratios were then averaged in consideration of the different levels of each parameter and are then plotted in Fig. 3.
In the Taguchi method, the S/N ratio gives a measure of the performance of a system or process by quantifying the relationship between the signal, which represents the desired output or target value, and the noise, which represents the variation from the target value.It measures how well the system meets the desired performance indicators (higherthe-better for SPF and GTC and lower-the-better for LCOH) in the presence of variability.A higher S/N ratio indicates a better quality of performance because the signal is more pronounced relative to the noise.Average S/N ratio plots (Fig. 3) help in selecting the best combination of parameter levels that lead to the desired outcomes with the least variability.The optimal parameter set ups for every response factor in all three configurations were chosen by looking into the level combinations that gave the highest average S/N ratios for every control factor (Fig. 3).
As shown in Fig. 3, the combination that gave the optimal SPF for configurations 1, 2, and 3 (parameter levels that give the highest S/N ratios) are A1B1C1D1E1F1G3H1, A3B1C1D1E2F1G3H1, and A3B1C1D1E3F1G3H1, respectively.The effect of a certain parameter on a performance indicator could be inferred from the difference between the highest and the lowest value of the S/N ratio.Here, it could be seen that the SPFs for all three configurations were notably more affected by the heat pump's discharge pressure (G) and the heat pump's outlet temperature (H).CO 2 heat pumps are known to exhibit an optimal discharge pressure that maximizes its efficiency [77].In the system considered here, the SPF-optimal high-side heat pump pressure seems to lie higher than 10 MPa.However, pressure levels higher than this were not explored anymore to avoid being too far from the operating conditions at which the CO 2 heat pump component model was calibrated.At a lower outlet temperature (H), the heat pump circulates a lower amount of CO 2 , reducing compressor work, and implements a larger change in specific enthalpy during the heat rejection process.These contributed to better SPFs.The parameter levels that optimized SPF for BHE spacing (B), BHE number (C), SC area (D), and SH return temperature (F) were seen to be consistent for all configurations, i.e., smaller BHE spacing (B), number (C), and SC area (D) and lower space heating return temperature (F) are preferrable for a better SPF for all configurations.However, note that making the BHE spacing (B), number (C), and SC area (D) small could entail a larger impact on the ground temperature, which could affect the long-term capability of the system to provide thermal energy.A smaller space heating return temperature (F) is expectedly more preferrable for all configurations since it implies a more efficient distribution system, which would require the circulation of less fluid to provide the SH demand.Applying the determined control factor combinations in simulation resulted in SPF values of 3.804, 4.025, and 4.018 for configurations 1, 2, and 3, respectively (Table 6).This shows that when the other performance indicators (LCOH and GTC) are disregarded, configuration 2 could give the best performance.Note that the SPFs here were calculated only after 1 year of operation.If the ground temperature changes significantly yearly, the annual SPF would change accordingly.
The combination that gave the optimal LCOH for configurations 1, 2, and 3 are A1B3C1D1E1F1G3H1, A1B2C1D1E2F2G3H1, and A1B2C1D1E3F1G3H1, respectively.The design parameters BHE length (A), BHE number (B), and SC area (C) were seen to be more influential to the LCOH compared to the other parameters.From this, it could be inferred that the capital cost contributes much to the LCOH.For all configurations, the LCOH is seen to be consistently insensitive to the BHE spacing (B).However, this is the result of the assumption in this study that the unit cost of the BHE is only related to its length and not the space it occupies.In reality, increasing the spacing may also increase costs, but possibly to a much lesser extent than when the BHE length is increased.Applying the determined control factor combinations in simulation resulted in LCOH values of 0.126 USD/kWh, 0.126 USD/ kWh, and 0.124 USD/kWh for configurations 1, 2, and 3, respectively.This shows that when the other performance indicators (SPF and GTC) are disregarded, the parallel configuration performs best (Table 6).

Table 3
The Taguchi L 27 (3 9 ) standard orthogonal array and the experimental plan.The combination that gave the optimal GTC for configurations 1, 2, and 3 are A3B3C3D3E1F3G1H3, A3B3C3D3E2F3G1H3, and A3B3C3D3E2F3G1H3, respectively.It could be seen that the GTC is most sensitive to the BHE (A, B, C) and SC (D) parameters.Higher values of these parameters could help in regulating the temperature reduction in the ground that could be induced by the system.However, increasing the BHE length (A), BHE number (C), and SC area (D) also increases the cost.Increasing the BHE spacing (B) might be a better way to manage ground thermal imbalance.Note though that the relationship between BHE spacing (B) and cost is artificially dissociated here.Applying the   determined control factor combinations in simulation resulted in GTC values of 100.257 %, 100.251 %, and 100.242 % for configurations 1, 2, and 3, respectively.This shows that when the other performance indicators (SPF and LCOH) are disregarded, the series configurations are preferable (Table 6) than the parallel.There are trade-offs when adjusting the values of the parameters in terms of their effects on the different performance indicators.Hence, it is necessary to perform multi-objective optimization to see the overall performances of the three configurations by just looking into one indicator (GRG) that cumulatively represents the three performance indicators.

Single-objective optimizations: Taguchi method -ANOVA
Performing ANOVA in Taguchi experiments is necessary for understanding the impact of different parameters and their levels on the variation in the output variables.It also helps in determining which factors are statistically significant and which are not.The results of the ANOVA are summarized in Table 7.When the P-value is less than the significance level (Significance codes: 0 '***', 0.001 '**', 0.01 '*', 0.05 '.'), it suggests that a parameter has a statistically significant impact on the S/N ratio and the performance indicator.The F value measures the ratio of the variance between groups to the variance within groups.A higher F value suggests a more significant impact of the factor.To calculate the contribution of a parameter to a performance indicator, the SS is considered.The SS measures the variability or dispersion of data points around a central value.It is used to quantify how much of the total variation in the output can be attributed to a certain parameter.
ANOVA was used to estimate the relative significance of each parameter in terms of percentage contribution to the different performance indicators (Table 7).The significance denotes which parameters could induce statistically significant effects on the response factor at different confidence levels.
The SPF was seen to be consistently most sensitive to the CO 2 heat pump's discharge pressure (G) for all three (3) configurations.The heat pump's outlet temperature (H) contributes significantly to the SPF too, especially for configurations 2 and 3.
The contribution of the BHE length (A) to the LCOH is rather substantial, given that the BHE is the most expensive component of the system.The BHE number (C) also contributed significantly but to a lesser extent since the levels chosen for this parameter were more constrained than the levels chosen for BHE length (A).The operation parameters heat pump's discharge pressure (G) and the heat pump's outlet temperature (H) played more significance to the LCOH of configurations 2 and 3 than it did for configuration 1.
The parameters BHE length (A), BHE spacing (B), BHE number (C), and SC area (D) consistently exhibited huge effects on the GTC of the three configurations since these are the parameters that can help disperse the thermal effects of the system to the ground.The SC area (D) was seen to be most significant to the GTC of configuration 1 while BHE spacing was seen slightly more important to configurations 2 and 3.

Multi-objective optimization: Taguchi-GRA
Taguchi-GRA was implemented to perform multi-objective optimization, with all performance indicators given equal weights.To perform this, normalization of the performance indicators was first done followed by the calculation of the GRGs.GRGs represent the overall performance of the system when multiple performance indicators are considered simultaneously.A value closer to 1 is preferred.The calculated GRGs, their respective S/N ratios, and rank are all given in Table 8.
The corresponding S/N ratios of the GRGs were calculated according to the "higher the better" concept, and their averages are plotted in Fig. 4. Taking note of the highest S/N ratios for each control parameter, the optimal combinations were noted to be A1B3C1D3E1F1G3H1, A1B3C1D1E2F1G3H1, and A1B3C1D1E3F1G3H1 for configurations 1, 2, and 3, respectively.The optimal parameter levels for parameters BHE length (A), BHE spacing (B), BHE number (C), SH return temperature (F), heat pump's discharge pressure (G), and heat pump's outlet temperature (H) are all identical for all the configurations considered.Keeping the BHE length (A), BHE number (B), SH return temperature (F), and outlet temperature (H) low but reasonable was seen as beneficial to the system.A shorter BHE length and using fewer BHEs lower capital expenses; low SH return temperature and heat pump's outlet temperature reduce the operations expenses (increase the SPF).A higher BHE spacing helps reduce the thermal impact of the system on the ground while operating near or at the optimal discharge pressure (G) increases the cost efficiency of the system.As for the area of the SC (D), it could be seen that a higher SC area is desirable for configuration 1, but smaller areas are preferred for configurations 2 and 3.It can also be inferred that consideration of low BHE length (A) and number (B) was only possible because of the heat provided by the SCs.
The desired SC-BHE circulation rate (E) varied in every configuration, with a smaller rate more beneficial to configuration 1, mid-level for configuration 2, and a higher rate for configuration 3.This shows the need to increase the circulation rate to facilitate enough heat absorption from the BHEs and SCs simultaneously when they are placed parallel to the heat pump, affecting the power consumption for pumping needed by the system.Looking at the magnitude differences between the highest and lowest S/N ratios for every control factor, it could be inferred that the BHE length (A) plays an important role in the overall performance of the system for all configurations, followed by the heat pump's discharge pressure (G) and outlet temperature (H).
Using the optimal parameter level combinations as inputs to simulations gave the optimized performance of every configuration (Table 9).The optimized configuration 1 gave the best overall performance among the three (GRG = 0.6875) while the optimized configurations 2 (GRG = 0.6644) and 3 (GRG = 0.6643) performed almost similarly.Nonetheless, looking into the GRGs obtained from the trial runs (Table 8), it could be seen that the parallel configuration produced more of the top 10 performers (GRG ≥ 0.60) in the orthogonal array.

Multi-objective optimization: Taguchi-GRA -ANOVA
ANOVA results (Table 10) show that the GRG of the 3 different configurations is most sensitive to the BHE length (A) and least sensitive to the SC area (D).This means that the BHE length should be kept as small as practical then the SC area can be increased to provide for the heating demand.The other relatively important parameters for all configurations include the CO 2 heat pump's discharge pressure (G) and the outlet temperature (H).It shows the significance of operating the heat pump at the optimal discharge pressure regardless of the   configuration of the system.The outlet temperature is set depending on the needs of the distribution system.The more efficient system will require a lower outlet temperature and, hence will contribute to a better overall performance.The SC-BHE mass flow (E) played a bigger role in configuration 1 than in the other set ups, while the BHE spacing (B) was seen as more important to the parallel configuration.

Conclusion
This study conducted an optimization analysis of the three different configurations (1: series A -SCs and BHEs in series with the fluid passing through the SCs first, 2: series B -SCs and BHEs in series with the fluid passing through the BHEs first, 3: parallel -SCs and BHEs in parallel with the fluid flow divided equally among them when solar energy is available) of a CO 2 SAGSHP system using the Taguchi method and GRA.Single-objective optimizations on the individual configurations were first implemented concerning the SPF, LCOH, and GTC using the Taguchi method.After that, multi-objective optimizations were performed using Taguchi-GRA considering the GRG, which represents the overall performance of the system.Here are the key findings and conclusions: • Multi-objective optimization showed that configuration 1 (Series A) could give the best overall performance when the SPF, LCOH, and GTC are considered simultaneously and given equal weights.In the system considered here, it gave a GRG of 0.6643, which translates to an SPF of 3.267, LCOH of 0.155 USD/kWh, and GTC of 100.021 %.These values demonstrate comparable SPF and LCOH relative to other heat pump systems, while the GTC slightly above 100 % suggests long-term thermal energy sustainability from the ground.• Optimization across different configurations shared some common parameter settings.These settings include keeping BHE length (A) short, maintaining high BHE spacing (B), using a low BHE number (C), keeping space heating return temperature (F) low, optimizing heat pump discharge pressure (G), and maintaining a low heat pump outlet temperature (H).Configuration 1 benefits from a high SC area, whereas Configurations 2 and 3 require lower SC areas.Additionally, the SC-BHE mass flow should be kept low for Configuration 1, somewhat higher for Configuration 2, and highest for Configuration 3.   • BHE length was found to be a crucial factor in overall performance, irrespective of the configuration.It is recommended to keep BHE length reasonably short and counterbalance the ground temperature decrease it would induce by increasing BHE spacing and SC area.• When optimizing solely for SPF, Configuration 2 delivered the best performance, achieving an SPF of 4.025.The SPFs for all three configurations were notably more affected by the heat pump's discharge pressure (G), which should be maintained at an optimal level, and the heat pump's outlet temperature (H), which should be kept as low as allowed.• Concerning only the LCOH, the parallel configuration outperformed the others, yielding an LCOH of 0.124 USD/kWh.The design parameters of BHE length (A), BHE number (B), and SC area (C) exerted a greater influence on LCOH compared to other parameters.Keeping these values low helps maintain a low LCOH.• Configuration 1 excelled when considering only GTC as a performance indicator to optimize, achieving a GTC of 100.267 %.GTC was most sensitive to BHE (A, B, C) and SC (D) parameters.Longer and more spaced BHEs, as well as larger SC areas, were shown to mitigate or reverse the thermal decline experienced by the ground during utilization.• Further optimization of system operations and comparisons with conventional subcritical HFC heat pumps under similar boundary conditions should also be explored.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The thermal demand used in this study and typical thermal demand in different seasons.

Fig. 3 .
Fig. 3.The average S/N ratios for each parameter and performance indicator.

Fig. 4 .
Fig. 4. The average S/N ratio for the GRGs for each configuration.

Table 1
Summary of parameters used for cost calculations.

Table 2
Parameters investigated and their levels.
*the flow rate of the fluid circulating through the SC, BHE, and the evaporator's water side.**thetemperature of the water as it comes out of the water side of the gas cooler.T.A.Sazon et al.

Table 4
Values of the performance indicators for the different system configurations.

Table 5 S
/N ratio.

Table 6
Results of the single-objective optimizations for each performance indicator in the different configurations.

Table 7
ANOVA tables for S/N ratios derived for each performance indicator for every configuration.

Table 8
Results of GRA.LCOH N GTC N SPF GRC LCOH GRC GTC GRC SPF N LCOH N GTC N SPF GRC LCOH GRC GTC GRC SPF N LCOH N GTC N SPF GRC LCOH GRC GTC GRC

Table 9
Results of the multi-objective optimization for every configuration.

Table 10
ANOVA for GRG-derived S/N ratios.