How to Evaluate the Operating Performance of Mid-Deep Geothermal Heat Pump Systems (MD-GHPs): A Study on a Multistage Evaluation Index System

Abstract

Mid-deep geothermal heat pump systems (MD-GHPs) use mid-deep borehole heat exchangers (MDBHEs) to extract heat from the geothermal energy at a depth of 2–3 km, and have been used for space heating in China over the last decade. This paper proposes a comprehensive and multilevel evaluation-index system to analyze and evaluate the energy performance of MD-GHPs. The multilevel evaluation index system consists of a target layer, a criterion layer, and an index layer, where the criterion layer is subdivided into six aspects and the index layer includes 26 specific indices, reflecting the geothermal resources, heat transfer performance of the MDBHEs, energy efficiency of the heat pump systems, building space heating demand, grid dynamic response capability, and energy-saving and economic benefits. Then, based on both expert survey results and case study data, the entropy weight method and the analytic hierarchy process are integrated to determine indicator weight coefficients among the multilevel evaluation indices, comprehensively considering both subjective and objective analyses. Furthermore, a fuzzy comprehensive evaluation model is conducted to integrate these weighted indices into a multi-criteria evaluation of MD-GHP performance. Finally, the proposed method was applied to evaluate the practical performance of four projects, returning scores of 61.56, 58.33, 72.73, and 78.41. These evaluations enable an overall assessment of the energy performance of MD-GHPs, reflecting the technical weaknesses and offering optimization guidance for system design and operation.

1. Introduction

Space heating for buildings is a major energy consumer in China. As of 2022, the heating area reached 16.7 billion m², consuming 217 million tons of standard coal and emitting 440 million tons of CO₂ [1]. Fossil fuels dominate, accounting for nearly 88% of heating sources, primarily through combined heat and power (CHP) generation, and gas or coal boilers, while electric boilers, heat pumps, and industrial waste heat recovery contribute only 12% [2]. This extensive reliance on fossil fuels leads to significant CO₂ emissions. On 22 September 2020, China announced its goal to peak CO₂ emissions by 2030 and achieve carbon neutrality by 2060. Consequently, adopting clean and low-carbon heating technologies is critical to meeting increasing demand [3].

Renewable energy-based heating technologies include solar thermal systems [4,5] and electric-driven heat pumps, such as ground source (GSHPs) and air source (ASHPs) heat pumps [6]. GSHPs are further classified into surface water heat pump systems (SWHPs), groundwater heat pump systems (GWHPs), and ground-coupled heat pump systems (GCHPs) [7]. GCHPs are considered more reliable and efficient as the ground temperature remains stable during the heating season, making them widely adopted in recent decades [8]. However, conventional GCHPs, with heat exchangers installed at depths below 200 m, face challenges including environmental and climatic sensitivities [9]. Moreover, large-scale land requirements [10] and thermal imbalance in the ground over time [11,12,13,14] limit their effectiveness. To address these limitations, GCHPs have been coupled with other renewable energy sources, such as ambient air [15,16], solar thermal energy [17,18,19], and biomass energy [20]. Nevertheless, these hybrid systems often increase land use and initial investment, restricting wider application.

A more direct solution involves mid-deep borehole heat exchangers (MDBHEs), which harness medium-depth geothermal energy from depths of 2 ~ 3 km, where temperatures range between 70 ~ 90 °C. Originally implemented in the U.S. and Europe, MDBHEs extract geothermal heat without withdrawing groundwater, as demonstrated by field tests (

Table 1. Field test results of the heat transfer performance of MDBHEs.

Area	Depth (m)	Inlet/Outlet Water Temp (°C)	Flow Rate (m3/h)	Heat Extraction Rate (kW)
Hawaii, USA [21]	1962	30/98	4.7	370
Penzlau, Germany [22,23]	2786	40/60	6.1	139
Weggis, Switzerland [24,25,26]	2300	32/40	10.7	100
Aachen, Germany [27]	2500	30/40	10.1	117

) [21,22,23,24,25,26,27].

Recently, this approach has been integrated with heat pump systems, termed mid-deep geothermal heat pump systems (MD-GHPs), for space heating in China. Field tests conducted by Deng revealed that for a single 2500 m DBHE with a ground-side water flow rate of 30 m³/h and inlet water temperature of 5 °C, the heat extraction rate could reach about 500 kW under continuous operation mode [28] and more than 600 kW under intermittent operation mode [29]. The heat pump’s coefficient of performance (COP) reached 5.43, while the overall system efficiency was 4.58. This high-temperature heat source enables MD-GHPs to outperform conventional shallow-depth systems in both energy savings and CO₂ reduction [28].

Research on MDBHEs has primarily focused on heat transfer performance, with numerical simulations identifying key influencing factors. Higher ground thermal conductivity and temperature, larger MDBHE depth and diameter, and larger thermal conductivity of the outer and lower thermal conductivity of the inner tubes enhance heat extraction [30,31]. Additionally, operational conditions, such as lower inlet water temperatures, higher flow rates, and intermittent operation could also optimize the heat extraction performance [32]. Long-term studies suggest that although ground temperature decreases over time, the reduction is less than 4.0% after 10 years, indicating stable long-term operation [33,34]. Deng [35] also studied the long-term performance of MDBHE arrays in single-line layout and revealed that with line space decreasing from 100 m to 10 m, the maximum accumulated heat extraction capacity per heating season decreased by 0.3% to 19.1% compared to a single MDBHE. Therefore, a line space larger than 30 m was recommended. Chen [36] compared deep enhanced U-tube borehole heat exchangers (EUBHEs) with coaxial tubes, finding that EUBHE systems achieve up to 1.2 MW of heat extraction per season, outperforming two DBHEs with equivalent borehole length. Previous studies have primarily focused on analyzing the heat transfer performance of MDBHEs, highlighting important optimization directions. However, due to the distinct operational conditions and performance characteristics of MD-GHPs compared to conventional shallow-depth geothermal heat pump systems, applying the same design parameters, equipment selection, control methods, or evaluation approaches used for traditional systems can lead to suboptimal energy performance [37]. As such, it is crucial to establish a comprehensive, multistage evaluation index system tailored to the unique attributes of MD-GHPs.

In the field of engineering, developing comprehensive evaluation systems has gained increasing attention. Researchers have applied various methods to construct comprehensive evaluation systems across various fields, supporting optimal decision-making and comprehensive benefit assessments. Among these methods, the analytic hierarchy process (AHP) method, proposed by Saaty and Vargas [38] in the mid-1970s, has been widely adopted due to its high feasibility and effectiveness, which is a systematic and hierarchical analysis method combining qualitative and quantitative analysis. For example, Yang et al. [39] developed an AHP-based evaluation system for ground source heat pumps and established grading standards for 15 evaluation indices. Meng et al. [40] combined the decision-making AHP with fuzzy mathematics to create a fuzzy comprehensive evaluation model for the green retrofitting of existing buildings. Man and Zhang [41] set up a green building evaluation system based on the AHP and grey clustering method. Yu et al. [42] developed an assessment method for green store buildings in China, which used the expert group decision AHP to determine the weight distributions of evaluation aspects highlighting the importance of indoor environmental quality, energy efficiency, and operation management within store buildings. Huang et al. [43] pointed out that in previous studies, an AHP was commonly used in planning fire station layouts considering critical factors affecting fire protection coverage. Additionally, they utilized a geographic information system to establish an optimal fire station layout. Ren et al. [44] established different operational models of the conventional geothermal heating system coupled with energy storage based on time-of-use electricity prices. They adopted a genetic algorithm to find the optimal decision variables minimizing the levelized cost of heat. Therefore, this study employed the AHP and entropy weight method (EWM), considering expert thinking and minimizing the influence of sample quality, to determine the weight distributions of evaluation indices. Linear reorganization was employed to integrate the subjective and objective weights. Additionally, fuzzy mathematics theory was used to decide each evaluation index rating.

As shown in Figure 1, this paper addresses the existing research gap by first introducing the field test methodology for MD-GHPs. Then, drawing upon a comprehensive review of the relevant literature, standards, and practical project experience, a hierarchical entropy weight method is applied to calculate the index weights. Subsequently, a fuzzy comprehensive evaluation model is constructed to integrate the results from multiple indices. A total of 26 evaluation indices are established, covering aspects such as the heat transfer performance of the MDBHEs, energy efficiency of the heat pump systems, terminal heat consumption, overall energy use, return on investment, carbon emissions, and other environmental impacts. The structural flow of the study is illustrated in Figure 1. This systematic evaluation framework aids in the optimization of MD-GHPs, from system design and equipment selection to practical operation. Additionally, by identifying weaknesses within the system, this evaluation method provides valuable insights for operational improvements and investment decisions, serving as a robust tool for enhancing the long-term sustainability and efficiency of MD-GHPs.

2. Establishment of Multistage Evaluation Index System of the MD-GHPs

Based on the on-site measurement, standard specification research [45,46,47], a literature review, and actual engineering surveys, this section develops a six-dimension and multi-level evaluation index system, which encompasses the following dimensions: Resource Conditions, Buried Pipe Heat Transfer Performance, Heating System Performance, Building Heating Needs, Grid Dynamic Response Capacity, and Energy Saving and Economic Benefits. The corresponding evaluation indices are identified to assess the application effectiveness of MD-GHPs. This comprehensive evaluation framework not only facilitates the assessment of existing engineering projects but also sets long-term operational performance targets for new installations, providing reference values for key indices. Ultimately, the framework serves as a guiding tool for system design, construction, and long-term operation.

2.1. Field Test Methodology of MD-GHPs

As shown in Figure 2, MDBHEs consist of an annulus tube with stainless-steel metal material and an inner tube with insulating material. During operation, the heat source water is pumped down through the annular tube, where it exchanges heat with the surrounding soil via the annulus tube wall. Then, it flows upward through the inner tube, where the heat is still transferred to annulus tube water as the inner tube wall is not completely adiabatic. Furthermore, the electric-driven heat pump system is integrated to raise the water supply temperature. Concurrently, the inlet water temperature of the MDBHE could be decreased to extract more geothermal energy.

Table 2. Basic information on DBHE and surrounding environment [35].

Parameters		Numerical Value
Annulus tube	External diameter (mm)	178
	Inner diameter (m)	160
	Material	Petroleum casing pipe (J55)
	Thermal conductivity coefficient (W/(m·°C))	54
	Specific heat capacity (J/(kg·°C))	0.47
	Density (kg/m3)	7820
Inner tube	External diameter (mm)	110
	Inner diameter (mm)	94
	Material	High-density polyethylene
	Thermal conductivity coefficient (W/(m·°C)	0.40
	Specific heat capacity (kJ/(kg·°C))	2.1
	Density (kg/m3)	930
Backfill material	Diameter (mm)	400
	Material	Grout
	Thermal conductivity coefficient (W/(m·°C)	2.0
	Specific heat capacity (kJ/(kg·°C))	0.85
	Density (kg/m3)	2700

lists the basic information on the MDBHE studied in this paper and its surrounding environment.

The heat transfer process in the MDBHE is analyzed using several key equations that describe the energy exchange between the system’s components:

(1) The heat transfer between the inner and annulus tube water is governed by the energy balance equation given as Equation (1),

$${\displaystyle \frac{\partial T_{w,i}}{\partial \tau }}+{\displaystyle \frac{\partial \left(u_{in}\times T_{w,i}\right)}{\partial z}}={\displaystyle \frac{K_i·\left(T_{w,o}-T_{w,i}\right)}{{\rho }_w·C_w{·A}_{in}}}$$

(1)

where $T_{w,i}$ and $T_{w,o}$ are the inner and annulus tube water temperatures, °C; $u_{in}$ is the inner tube water velocity, m/s; $A_{in}$ is the inner tube flow area, m²; $\tau$ is the time, s; z is depth parameter, m; and $K_i$ is the inner tube thermal conductivity per unit length, W/(m·K), calculated using Equation (2),

$$K_i=\pi /({\displaystyle \frac{1}{h_1·d_{ii}}}+{\displaystyle \frac{1}{2{·\lambda }_{in}}}\ln {\displaystyle \frac{d_{io}}{d_{ii}}}+{\displaystyle \frac{1}{h_2·d_{io}}})$$

(2)

where $d_{ii}$ and $d_{io}$ are the inner and outer diameters of the inner tube, m; ${\lambda }_{in}$ is inner tube thermal conductivity, W/(m·K); $h_1$ is the inner tube inner surface heat transfer coefficients; and $h_2$ is the inner tube outer surface heat transfer coefficients, W/(m²·K).

(2) The heat transfer between the annulus tube water and the ground borehole wall is described using Equation (3),

$${\displaystyle \frac{\partial T_{w,o}}{\partial \tau }}+{\displaystyle \frac{\partial \left(u_o{·T}_{w,o}\right)}{\partial z}}={\displaystyle \frac{K_i·(T_{w,i}-T_{w,o})}{{\rho }_w·C_w·A_o}}+{\displaystyle \frac{K_o·(T_b-T_{w,o})}{{\rho }_w{·C}_w·A_o}}$$

(3)

where $u_o$ is annulus tube water velocity, m/s; T_b is the ground borehole wall temperature, °C; and $K_o$ is annulus tube thermal conductivity per length, W/(m·K), calculated using Equation (4),

$$K_o=\pi /({\displaystyle \frac{1}{h_3·d_{oi}}}+{\displaystyle \frac{1}{2·{\lambda }_o}}\ln {\displaystyle \frac{d_{oo}}{d_{oi}}}+{\displaystyle \frac{1}{2·{\lambda }_g}}\ln {\displaystyle \frac{d_b}{d_{oo}}}+R_c)$$

(4)

where $d_{oo}$ and $d_b$ are the diameters of the annulus tube and the borehole wall, m; ${\lambda }_o$ is annulus tube thermal conductivity, W/(m·°C); $h_3$ is the inner surface heat transfer coefficient of the annulus tube, W/(m²·K); and $R_c$ represents the contact thermal resistance between MDBHE and backfill material, (m·K)/W.

(3) The heat transfer in the surrounding ground is expressed using Equation (5),

$${\rho }_g·C_g·{\displaystyle \frac{\partial T_g}{\partial \tau }}={\displaystyle \frac{1}{r}}·{\displaystyle \frac{\partial }{\partial r}}\left(r{·\lambda }_g·{\displaystyle \frac{\partial T_g}{\partial r}}\right)+{\displaystyle \frac{\partial }{\partial z}}({\lambda }_g·{\displaystyle \frac{\partial T_g}{\partial z}})$$

(5)

where T_g is ground temperature, °C; ${\rho }_g$ is ground density, kg/m³; $C_g$ is ground heat capacity, kJ/(kg·K), and ${\lambda }_g$ is ground thermal conductivity, W/(m·K).

These equations are used to analyze the heat exchange processes in the MDBHE system. They provide a detailed understanding of the thermodynamic behavior of the system, including the interaction between the water, the tubes, and the surrounding ground.

The practical energy performance of MD-GHPs in four projects has been measured and analyzed.

Table 3. Basic Information on Field Test Systems.

Project	Building Function	Space Heating Area (m2)	Rated Heating Capacity (kW)	Indoor Terminals	Depth of MDBHE (m)	Numb of MDBHE	Monitoring Period
MG-1	Residence	56,000	2600	Radiant floor	2000	4	2 heating seasons
MG-2	Residence	133,400	5680	Radiant floor	2500	8	2 heating seasons
MG-3	Residence	185,100	7560	Radiant floor	2500	10	2 heating seasons
MG-4	Office	33,160	2410	FCU	2800	2	2 heating seasons

outlines the basic information of these projects. Three projects were residential buildings, while one was an office building. The heating areas ranged from 30,000 to 180,000 square meters. The residential projects utilized floor radiant heating, whereas the office building employed fan coil units. Based on the actual heating demand of the projects, 2 to 10 MDBHEs were designed, each with depths exceeding 2000 m.

Figure 3 shows the measuring points used during the field tests of MD-GHP projects. Since 2015, the research team has conducted field tests and continuous monitoring in multiple practical projects to assess the operational performance of MD-GHPs. Following the test methods and guidelines outlined in the standard [48], key monitoring data includes user-side circulation flow, supply and return water temperatures, water system pressure distribution, heat-source-side circulation flow, inlet and outlet water temperatures, and power consumption of major equipment such as heat pumps and water pumps. To ensure accuracy, sensors with appropriate precision and range were selected based on error analysis; the instruments used are detailed in

Table 4. Field test instrument information.

Physical Quantities	Measuring Accuracy	Measuring Range
Temperature sensor	±0.1 °C	−50~100 °C
Clamp-on ultrasonic flow meter	±1%	0~30 m/s
Three-phase electric power meter	±0.5%	0~800 kW

. In addition, regular single-point testing was performed alongside long-term monitoring to cross-validate the data, reducing random errors and enhancing the reliability of temperature readings.

2.2. Construction of the Multistage Evaluation Index System

The multistage evaluation index system comprises three hierarchical levels: the target layer, criterion layer, and index layer. The target layer defines the primary objectives and overarching purpose of the MD-GHPs. The criterion layer is further subdivided into six main aspects, which reflect the geothermal resource availability, heat transfer performance of the MDBHEs, and energy efficiency of the heat pump systems, as well as economic and energy-saving impacts.

The index layer contains 26 specific evaluation indices, each measured through varying approaches. Some indices are determined by consulting relevant literature and location-based data. Others are derived through short-term or long-term measurements conducted on-site. Additionally, certain indices require calculations based on measured data. Among these indices, some parameters, like the heating performance coefficient, directly influence the system or MDBHE performance, while others evaluate user-side energy savings or economic outcomes, indirectly reflecting system effectiveness. Figure 4 illustrates the hierarchical structure of this evaluation index system.

2.2.1. Criterion Layer-A: Resource Conditions

The criterion layer-A focuses on the key meteorological and geological resources, which could be classified into the following parts:

• A1: Ground Surface Temperature (T_gsur), which is a key indicator for evaluating ground source heat pump systems, influenced by solar radiation and ambient temperature. This paper adopts a constant boundary condition, with the temperature set as the annual average. Evaluation standards are based on the 2023 data distribution for major Chinese cities [49].
• A2: Soil Geothermal Gradient (T_gd), which reflects the temperature increase per 100 m below the Earth’s surface. This is linked to the Earth’s geothermal energy, with deeper soil layers (beyond 200 m) exhibiting an average gradient of 3 °C per 100 m [50]. This gradient directly influences the heat extraction efficiency of MDBHEs. Figure 5 outlines the average soil geothermal gradients across China’s climate regions, forming the basis of the evaluation standard [51]. Regions with low heating demand, such as those with hot summers and warm winters, are excluded from the study.

3.

A3: Soil Thermal Conductivity (λ_g), which quantifies the soil’s ability to conduct heat, and is influenced by the composition of dry soil, gas, and pore fillers [52]. Measurement methods are broadly categorized as steady-state [53,54] and transient-state [55,56,57].

Table 5. Ground thermal conductivity values of typical rocks and soils.

Soil and Rock		λg (W/(m·°C))
Soil	Gravel	0.98
	Coarse Sand	1.02
	Fine Sand	1.03
	Silty Sand	2.07
	Clay	1.25
Rock	Limestone, Dolomite	2.46
	Water-Eroded Limestone	3.56
	Sandstone	4.5
	Shale	2.53
	Cracked Igneous and Metamorphic Rocks	4.61
	Uncracked Igneous and Metamorphic Rocks	4.59

summarizes typical ground thermal conductivity, forming the basis for this indicator’s evaluation standard [58,59].

4.

A4: Soil Density (ρ_g) is also a crucial parameter. This paper compiles density data from the standards presented in

Table 6. Typical density values of rocks and soils.

Rock	ρg (kg/m3)
Granite	2700
Limestone	2700
Sandstone	2600
Calcium Sand (Moisture Content of 43%)	1670
Dry Quartz Sand (Medium-Fine Grained)	1650
Quartz Sand (Moisture Content of 8.3%)	1750
Sandy Clay (Moisture Content of 15%)	1780

[58], establishing evaluation standards based on these data for further analysis.

5.

A5: Soil Heat Capacity (C_g), which is influenced by moisture content, density, and texture. The specific heat capacities of sandy soils can be assessed through sandbox experiments [52]. This paper uses the thermal capacities of common underground rocks and soils presented in

Table 7. Typical heat capacity values of rocks and soils.

Rock	Cg (J/(kg °C))
Granite	794
Limestone	920
Sandstone	878
Calcium Sand (Moisture Content of 43%)	2215
Dry Quartz Sand (Medium-Fine Grained)	794
Quartz Sand (Moisture Content of 8.3%)	1003
Sandy Clay (Moisture Content of 15%)	1379

, along with evaluation standards for these measurements [58].

6.

A6: Underground Seepage (P_e). Groundwater exists in various forms [60], which can significantly influence the thermal and physical properties of soil, thereby impacting the heat transfer performance of MDBHEs. This paper introduces the Peclet Number (P_e), a dimensionless parameter that quantifies the ratio of convection rate to diffusion rate in the context of underground seepage. It compares the intensity of thermal convection to thermal conduction and is expressed mathematically as Equation (6),

$$P_e={\displaystyle \frac{{\rho }_{gw}{C}_{gw}{u}_{gw}L}{{\lambda }_{eff}}}$$

(6)

where ${\rho }_{gw}$ is the groundwater density, kg/m³, $C_{gw}$ is the groundwater heat capacity, J/(kg·K); $u_{gw}$ is the groundwater flow rate through unit cross-sectional area, m/s·m²; $L$ refers to a fixed length, generally related to the length of the pipe, m; and ${\lambda }_{eff}$ is the effective thermal conductivity, W/(m·K).

The Peclet Number (P_e) varies significantly in underground rock and soil, with reference values provided for different geological conditions in

Table 8. Typical underground seepage Peclet Number of rocks and soils.

Soil and Rock		Pe
Soil	Gravel	5.72 × 102
	Coarse Sand	1.34 × 10
	Fine Sand	1.15
	Silty Sand	1.28 × 10−2
	Clay	3.24 × 10−5
Rock	Limestone, Dolomite	5.29 × 10−3
	Water Eroded Limestone	5.28
	Sandstone	1.77 × 10−3
	Shale	1.05 × 10−6
	Cracked Igneous and Metamorphic Rocks	6.32 × 10−2

[58,59]. This study also presents evaluation criteria for underground seepage to analyze the implications of groundwater movement.

2.2.2. Criterion Layer-B: Heat Transfer Performance of MDBHEs

This criterion layer conducts an in-depth analysis of the main factors affecting the heat transfer performance of MDBHEs. By exploring various influences on heat transfer rates, it proposes essential evaluation indices to optimize design and enhance overall heat transfer performance.

• B1: Transient Heat Extraction Power (Q_e) per MDBHE is a crucial indicator for evaluating system performance, calculated using Equation (7) [61],

  $$Q_e=\left(17.61\times T_{gd}\times {\lambda }_g+49.2\times T_{gd}-2.23\times {\lambda }_g\right)\times {\displaystyle \frac{H}{1000}}-8.63\times T_{gd}\times {\lambda }_g-61.93\times T_{gd}-7.92$$

  (7)

  where $H$ denotes the installation depth of the heat exchanger, m.
• B2: Comprehensive Heat Transfer Resistance (KF). Achieving higher outlet water temperatures and greater heat extraction power is essential for optimal MDBHE performance. However, these objectives can be inversely related—lowering inlet temperatures while increasing circulating water flow rate may enhance heat extraction but can also reduce outlet temperatures. To evaluate heat extraction performance, this study introduces the concept of comprehensive heat transfer resistance (KF) based on Entransy theory. It is calculated using Equation (8),

  $$KF={\displaystyle \frac{\mathsf{\Delta }{E}_n}{Q^2}}={\displaystyle \frac{\mathsf{\Delta }T}{Q}}$$

  (8)

  where $\mathsf{\Delta }T$ represents the temperature difference, °C, obtained from actual inlet and outlet temperatures, and $Q$ is the heat extraction power, kW. A lower KF indicates improved heat exchange with the surrounding soil; however, this must be balanced to avoid excessive heat dissipation from the inner pipe water. In practical scenarios, where measuring internal pipe water temperature distribution is challenging, using actual inlet and outlet temperatures provides a straightforward means of evaluating how heat transfer resistance impacts the ground heat exchanger’s performance.
• B3: Cumulative Heat Extraction Capacity (Q_long), which represents the total heat extracted per MDBHE over a total heating season. This serves as a key indicator for assessing the system’s long-term heat extraction capability and operational efficiency. The cumulative value is typically obtained through continuous monitoring over a prolonged period.
• B4: Long-term Operating Performance (OP_long). This paper uses annual temperature attenuation (°C/year) to express long-term performance. Overuse of an MDBHE can cause a marked ground temperature decrease, leading to a reduction in the heat transfer performance of the MDBHE. By continuously monitoring the average outlet water temperature of the underground pipes, this indicator can be effectively tracked. Excessive temperature attenuation directly impacts the long-term viability of MD-GHPs. Hence, maintaining an optimal balance between heat extraction and temperature stability is critical for ensuring sustained system performance. The evaluation criteria for this indicator are based on practical measurements from several projects and a comprehensive literature analysis.

2.2.3. Criterion Layer-C: Energy Performance of MD-GHPs

MD-GHPs consist of three main components: the heat pump unit, the user-side water system, and the heat-source-side water system. The performance of the heating system is directly influenced by the efficiency of each component. This criterion layer selects the performance coefficient of the heat pump unit and the transport performance coefficients of both the user-side and heat source-side water systems as key indices for analysis.

• C1: Heating Performance Coefficient of Heat Pump (COP_hp). The energy efficiency of the heat pump unit is crucial for MD-GHPs. The unit heating performance coefficient (COP_hp) is a dimensionless parameter representing the ratio of heating capacity to power consumption, thus reflecting the energy efficiency of the heat pump unit. Additionally, the system performance coefficient (COP_sys) expands this evaluation by incorporating the water systems. The relationships are expressed as Equations (9)~(11),

  $${\mathrm{C}\mathrm{O}\mathrm{P}}_{hp}={\displaystyle \frac{Q_s}{N_{hp}}}$$

  (9)

  $${\mathrm{C}\mathrm{O}\mathrm{P}}_{sys}={\displaystyle \frac{Q_s+Q_{h2}-Q^{\prime }}{N_s}}$$

  (10)

  $$N_s=N_{hp}+N_u+N_h$$

  (11)

  where $Q_s$ is the heating capacity of the MD-GHPs, kW; $Q_{h2}$ is the heating value of the user-side water pump, which comes from the electricity input and can generally be ignored, kW; $Q^{\prime }$ is the heat dissipation of pipelines and equipment, which can also generally be ignored, kW; and $N_s$ is the total input power of MD-GHPs, including the heat pump unit ($N_{hp}$), the user-side water pumps ($N_u$), and the heat-source-side water pumps ($N_h$), kW. Evaluation standards for COPhp are based on established regulations [61], as well as empirical data from engineering experience.
• C2: User-side water system transport performance coefficient (WTF_u) is a dimensionless metric that quantifies the efficiency of heat delivery, and can be calculated using Equation (12),

  $${\mathrm{W}\mathrm{T}\mathrm{F}}_{user}={\displaystyle \frac{Q_s+Q_{h2}-Q^{\prime }}{N_u}}$$

  (12)
• C3: Heat-source-side water system transport performance coefficient (WTF_h) is calculated similarly using Equation (13),

  $${\mathrm{W}\mathrm{T}\mathrm{F}}_{heat}={\displaystyle \frac{Q_s}{N_h}}$$

  (13)

2.2.4. Criterion Layer-D: Building Space Heating Demand

The evaluation of space heating demand in building systems is crucial for optimizing energy efficiency and ensuring occupant comfort. Understanding the specific demands on heating systems allows for better alignment of resources and technologies, directly influencing system design and performance. This layer focuses on key indices that collectively assess the heating demand of buildings, incorporating factors such as heating excess rate, transmission heat loss, and thermal storage performance.

• D1: Heating Demand (HD) quantifies the heat load required for a building, calculated by multiplying the building’s heating heat index by its area over a heating period or unit time, W/m². This indicator serves as a foundation for understanding the thermal requirements of a building and is often based on established thermal index standards that account for energy-saving measures. Evaluation standards are developed from practical engineering experience [61].
• D2: Heating Excess Rate (HER) is a dimensionless quantity that reflects inefficiencies in heating distribution caused by inadequate user-side regulation [61]. This value can be determined through detailed field tests, simulations, or empirical assessments based on the scale of the heating area. By correlating heating demand with excess heating, this indicator helps identify areas for improvement in energy distribution.
• D3: Transmission Heat Loss (THL), which is also dimensionless, indicates the heat lost during distribution within the heating network [61]. Measured by dividing the system’s heat loss by the distributed heat supply, this indicator highlights the importance of minimizing losses in achieving efficient heating performance. When direct measurement is not feasible, THL can be estimated through testing or calculations from heating station data.
• D4: Thermal Storage Performance (TSP) assesses the building’s ability to maintain indoor temperatures within an acceptable range while considering its heat storage characteristics. This involves defining an allowable temperature range (e.g., 18–24 °C in winter) and the duration over which this range can be sustained (e.g., for a minimum of 5 h). This indicator reflects the building’s capacity to absorb and retain heat, enhancing overall comfort and energy efficiency. TSP values can be obtained through actual testing, allowing for empirical validation of thermal performance.

2.2.5. Criterion Layer-E: Grid Dynamic Response Capacity

This criterion layer evaluates the system’s capability to respond dynamically to power grid demands, focusing on equipment efficiency, clean energy integration, and overall system resilience. Key indices reflect how well the system can utilize clean energy sources and adapt to operational conditions while maintaining economic viability.

• E1: Equipment Utilization Rate (UR_equ) quantifies the ratio of the peak load achieved during a full year of operation of an MD-GHP project to the system’s designed installed capacity. It is defined using Equation (14),

  $${\mathrm{U}\mathrm{R}}_{\mathrm{e}\mathrm{q}\mathrm{u}}={\displaystyle \frac{{HL}_{act-50}}{{HL}_{design}}}$$

  (14)

  where ${HL}_{act-50}$ refers to the peak load value during the full-year operation without guaranteeing 50 h, kW; and ${HL}_{design}$ refers to the value of the system’s designed installed capacity, kW. A UR_equ close to 100% signifies optimal utilization of the system’s capacity, minimizing waste from over-design. Evaluation criteria for UR_equ are based on practical testing and a literature review, promoting effective capacity utilization in actual projects.
• E2: Clean Energy Utilization Rate (EUR_Clean) measures the proportion of clean energy utilized in the total energy consumption of an MD-GHP project over a specified time frame. This indicator assesses the environmental sustainability of the energy usage structure by reflecting the extent to which clean energy contributes to overall energy consumption [62]. It is defined using Equation (15),

  $${\mathrm{E}\mathrm{U}\mathrm{R}}_{\mathrm{C}\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{n}}={\displaystyle \frac{E_{Clean}}{E_{total}}}$$

  (15)

  where $E_{Clean}$ refers to the total amount of clean energy used during system operation, which can be in the form of heat or electricity, kW; and $E_{total}$ refers to the total energy consumption during system operation, which can be in the form of heat or electricity, kW. EUR_Clean provides insights into the project’s environmental impact, allowing for a direct evaluation of clean energy integration and its contribution to sustainability goals.
• E3: Cheap Energy Utilization rate (EUR_Cheap) quantifies the proportion of energy costs incurred under relatively cheaper charging conditions in relation to the total energy costs during operation over a defined time scale. This indicator highlights the economic advantages of utilizing cost-effective energy sources. It is defined using Equations (16) and (17),

  $$OC=\sum _{i=1}^n\left(E_i{EP}_i\right)$$

  (16)

  $${\mathrm{E}\mathrm{U}\mathrm{R}}_{\mathrm{C}\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{p}}={\displaystyle \frac{{OC}_{Cheap}}{OC}}$$

  (17)

  where $OC$ refers to the total energy consumption costs generated by buildings, CNY. $E_i$ refers to the consumption of the i-th type of energy in the building, with units determined by the type of energy; ${EP}_i$ refers to the unit price of the i-th type of energy, which can be determined based on local market prices; and ${OC}_{Cheap}$ refers to the energy costs generated during system operation under relatively cheaper charging conditions, which can be directly obtained through energy metering, CNY. EUR_Cheap reflects the cost-effectiveness of the energy used, demonstrating the financial sustainability of the system while highlighting opportunities for optimizing energy expenditures.
• E4: System Resilience Coefficient (SRC) serves as an indicator to evaluate the system’s capacity for peak shaving within the power grid. It is defined as the effective duration of peak shaving participation in the power grid. For instance, during periods of low electricity demand or surplus renewable energy electricity, energy is stored within the building itself or with energy storage capabilities. Conversely, during periods of high electricity demand or a shortage of renewable energy electricity, energy is discharged from the building or energy storage capabilities to effectively reduce electricity load, thereby actively participating in peak shaving for the power grid. It is defined using Equation (18),

  $$SRC={\sum }_{t=t_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}1}}^{t_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}2}}\left(t|N_{\mathrm{c}\mathrm{u}\mathrm{r},\mathrm{t}}\ge 90\%N_{\mathrm{s}\mathrm{u}\mathrm{b}}\right)\cdot ∆t$$

  (18)

  where $t_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}1}$ is the starting time of the power grid peak shaving demand response event; $t_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}2}$ is the ending time of the power grid peak shaving demand response event; $N_{\mathrm{c}\mathrm{u}\mathrm{r},\mathrm{t}}$ is the peak load reduction of the heat pump system during the peak shaving response period, kW; $N_{\mathrm{s}\mathrm{u}\mathrm{b}}$ is the load reduction target value during the power grid peak shaving demand response event, kW; and $∆t$ is the time step, typically set as 1 h, but can also be adjusted to different intervals such as 0.5 h, 0.25 h, etc.

2.2.6. Criterion Layer-F: Energy Saving and Economic Benefits

Criterion Layer-F focuses on the energy-saving benefits facilitated by MD-GHPs and the accompanying economic advantages. This layer integrates both environmental and economic indices as follows:

• F1: Operational CO₂ Emission (CE_s) quantifies the emissions generated by an MD-GHP per unit area. The calculation, based on the consumption of various energy types and their respective carbon emission factors, is governed by relevant standards [62]. This paper further refines the evaluation criteria using empirical data, ensuring the standards align with actual performance metrics using Equations (19) and (20),

  $$CE=\sum _{i=1}^n\left(E_i{EF}_i\right)$$

  (19)

  $${CE}_s={\displaystyle \frac{CE}{A}}$$

  (20)

  where $CE$ refers to the total amount of carbon dioxide emissions generated during the operation of an MD-GHP, kg; ${EF}_i$ refers to the carbon emission factor of the i-th type of energy, which can be determined based on relevant standards; and $A$ refers to the space heating area, m².
• F2: Operational Primary Energy Cost (E_s,pri) reflects the total primary energy consumed during the operation of MD-GHPs per unit area. Primary energy refers to energy derived directly from natural sources without processing, with renewable sources typically considered. By establishing evaluation criteria based on area-based metrics, this indicator facilitates comparisons across various buildings, enhancing understanding of primary energy utilization efficiency. It is defined using Equation (21),

  $$E_{s,pri}={\displaystyle \frac{E_{pri}}{A}}$$

  (21)

  where $E_{pri}$ refers to the total amount of primary energy consumed in the system using MD-GHPs, which can be directly obtained through measurement, with units determined by the type of energy.
• F3: System energy saving rate (ESR) measures energy savings achieved by MD-GHPs compared to traditional energy systems. It reflects the difference in energy consumption under identical load conditions, allowing for a comprehensive analysis of energy-saving effects. The study formulates evaluation criteria based on literature and practical project data, which aid in assessing the effectiveness of energy-saving measures. It is defined using Equation (22),

  $$ESR={\displaystyle \frac{E_{other}-E_{total}}{E_{other}}}$$

  (22)

  where $E_{other}$ refers to the total energy consumption generated by other types of systems under the same load conditions, kgce; and $E_{total}$ refers to the total energy consumption generated by MD-GHPs, kgce.
• F4: System operating cost (OC_s) assesses the total operational costs of MD-GHPs per unit area, providing valuable insights into the cost-effectiveness of the system. By comparing these costs across different projects, stakeholders can better understand the financial implications of adopting MD-GHPs. Evaluation criteria are derived from practical project data and the relevant literature, ensuring a robust framework for analysis [63,64]. It is defined by Equation (23),

  $${OC}_s={\displaystyle \frac{OC}{A}}$$

  (23)
• F5: Static investment payback period (SIPP) is calculated by dividing the incremental investment cost of implementing MD-GHPs by the annual operational cost savings. A payback period of 10 years is considered reasonable, providing a benchmark for evaluating the return on investment. It is defined with Equation (24),

  $$\mathrm{S}\mathrm{I}\mathrm{P}\mathrm{P}={\displaystyle \frac{C_{\mathrm{M}\mathrm{D}-\mathrm{G}\mathrm{H}\mathrm{P}\mathrm{s}}}{{OC}_{\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}}-{OC}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}}$$

  (24)

  where $C_{\mathrm{M}\mathrm{D}-\mathrm{G}\mathrm{H}\mathrm{P}\mathrm{s}}$ refers to the incremental investment cost of MD-GHPs, CNY; ${OC}_{\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}}$ refers to the total operational energy consumption costs of a gas boiler system, CNY; and ${OC}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$ refers to the total operational energy consumption costs of MD-GHPs, CNY.

3. Analytic Hierarchy Process–Fuzzy Comprehensive Evaluation

Selecting an appropriate weighting method for evaluating comprehensive benefits in engineering is critical. The analytic hierarchy process (AHP) is a multi-criteria decision-making technique that decomposes complex decisions into smaller, manageable components. AHP promotes expert consistency through rigorous testing, resulting in reasonable weighted outcomes. Conversely, the entropy weight method (EWM) assigns weights based on the inherent information content of criteria, minimizing subjective bias [65,66,67,68]. In this section, the EWM and the AHP are integrated based on the constructed index system, which comprehensively considers both subjective and objective analyses. Then, a fuzzy comprehensive evaluation model, along with a multi-level evaluation index system tailored for MD-GHPs, is established through the fuzzy relationship matrix and evaluation factor, as illustrated in Figure 6.

3.1. Construction of the Entropy Weight Method

EWM provides a systematic approach to weight allocation when faced with numerous criteria, assigning weights based on the entropy of each criterion’s information. Entropy serves as a measure of uncertainty, and higher entropy indicates greater uncertainty. Consequently, weights are assigned inversely to entropy values, with lower entropy corresponding to higher importance. The effectiveness of EWM in facilitating objective decision-making is especially pronounced when prioritizing criteria based on their significance. The steps in the MD-GHP Multistage Evaluation Index System are as follows:

• Constructing a Standardized Matrix: The MD-GHP data matrix, comprising n evaluation indices and m evaluation projects, is normalized to form a standardized matrix using Equations (25)–(27),

  $$Y={\left(y_{ij}\right)}_{n\times m}$$

  (25)
  Then, for indices where smaller values are preferable:

  $$y_{ij}={\displaystyle \frac{max\left(x_{ij}\right)-x_{ij}}{max\left(x_{ij}\right)-min\left(x_{ij}\right)}}$$

  (26)
  While, for indices where larger values are preferable:

  $$y_{ij}={\displaystyle \frac{x_{ij}-min\left(x_{ij}\right)}{max\left(x_{ij}\right)-min\left(x_{ij}\right)}}$$

  (27)

  where $0\le y_{ij}\le 1$; $x_{ij}$ refers to the value of the i-th indicator of the j-th evaluation project; $max\left(x_{ij}\right)$ is the maximum value in $x_{ij}$; and $min\left(x_{ij}\right)$ is the minimum value in $x_{ij}$.
• Entropy Calculation: The entropy value H_i for the i-th indicator is calculated using Equations (28)–(30),

  $$H_i=-k\sum _{j=1}^m{f}_{ij}ln{f}_{ij}\left(i=\mathrm{1,2},\cdots ,n;j=\mathrm{1,2},\cdots ,m\right)$$

  (28)

  $$f_{ij}={\displaystyle \frac{y_{ij}}{\sum _{j=1}^m{f}_{ij}}}$$

  (29)

  $$k={\displaystyle \frac{1}{lnm}}$$

  (30)
  Additionally, to prevent situations where $f_{ij}=0$, a correction formula is introduced as shown in Equation (31),

  $$f_{ij}={\displaystyle \frac{{1+y}_{ij}}{\sum _{j=1}^m{(1+f}_{ij})}}$$

  (31)
• Weight Calculation: after determining the entropy values, weights are computed using Equation (32),

  $$w_i={\displaystyle \frac{1-H_i}{n-\sum _{i=1}^n{H}_i}}$$

  (32)

  where $\sum _{i=1}^n{w}_i=1$, $0\le w_i\le 1$, and $w_i$ refers to the weight value of each evaluation indicator.

3.2. Construction of the Analytic Hierarchy Process

AHP is a decision-making methodology that integrates qualitative and quantitative approaches to address complex problems. It structures decision criteria hierarchically, using pairwise comparisons to form a judgment matrix, which determines priority weights through mathematical algorithms. These weights reflect a criterion’s relative importance in decision objectives. AHP emphasizes consistency evaluation to ensure reliable judgments and adjust for any inconsistencies. This systematic approach enhances clarity and confidence in navigating complex decisions.

• Judgment Matrix Construction: This paper employs the 1-9 scale method for pairwise comparisons, quantifying the relative importance of each indicator in a hierarchical structure. The scale includes five basic categories: Equal (1), Weak (3), Moderate (5), Strong (7), and Absolute (9), with intermediate values (2, 4, 6, 8) providing finer gradations between them.

  Table 9. Evaluation of 1-9 scale method.

  Evaluation Scale	Definition	Explanation
  1	Equal Strong	Both are equally important
  3	Weak Strong	Former is weakly important than latter
  5	Strong	Former is obviously more important than latter
  7	Very Strong	Former is more important than latter
  9	Absolution	Former is absolutely more important than latter
  2, 4, 6, 8	Intermediate	Intermediate value between the above levels

  illustrates the evaluation scale. Based on this scale, a judgment matrix A can be constructed for further analysis.

2.

Weight Analysis: Once the judgment matrix A is established, the m-th root of each row’s product is computed to form an m-dimensional vector. The vector is then normalized to determine the weights of each indicator. The maximum eigenvalue λ_max is calculated using the weight matrix via Equations (33)–(35),

$${\dot{w}}_i=\sqrt[m]{\prod _{j=1}^m{a}_{ij}}$$

(33)

$$w_i={\displaystyle \frac{{\dot{w}}_i}{\sum _{j=1}^m{\dot{w}}_j}}$$

(34)

$${\lambda }_{max}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\left(Aw\right)}_i}{w_i}}$$

(35)

where ${\dot{w}}_i$ refers to the weight without normalization; and $w_i$ refers to the normalized weights.

3.

Consistency Calculation: to ensure the validity of the assigned weights, the consistency index (CI) and consistency ratio (CR) are calculated using Equations (36) and (37),

$$CI={\displaystyle \frac{{\lambda }_{max}-n}{n-1}}$$

(36)

$$CR={\displaystyle \frac{CI}{RI}}$$

(37)

where $CI$ refers to the consistency index; $RI$ refers to a random index; $CR$ refers to the consistency ratio; and $n$ refers to the order of the judgment matrix. If $CR<0.1$, this indicates that the constructed judgment matrix has passed the consistency check, and the assignment of weights to indices is reasonable; otherwise, if $CR>0.1$, this indicates that the consistency check has not been passed, and it is necessary to readjust the judgment matrix and recalculate the weights.

3.3. Integration of the Hierarchical and Entropy Weights

Using either the EWM or the AHP independently can result in either overly objective or overly subjective weights. To address this limitation, this paper combines both methods to calculate the comprehensive weight, $w_{com}$, using Equation (38) [66],

$$w_{com}=a\times w_{sub}+\left(1-a\right)\times w_{obj}$$

(38)

where $w_{sub}$ refers to the objective weight obtained from the EWM; $w_{obj}$ refers to the subjective weight obtained from the AHP; and $a$ refers to the combination coefficient of the two weights, with a range from 0 to 1. Here, it is set as 0.5, indicating that both weights are considered equally important.

3.4. Construction of a Fuzzy Comprehensive Evaluation Model

A comprehensive fuzzy evaluation model is essential for assessing the feasibility and performance of MD-GHPs. This model integrates both objective and subjective inputs, offering a systematic framework for decision-making. When combining the analytic hierarchy process and entropy weight method, AHP assigns priority weights through pairwise comparisons, while EWM objectively allocates weights based on criteria entropy, capturing uncertainty. The combination of these methods addresses the complexity of evaluation and improves both analytical depth and practical relevance. Fuzzy logic principles account for uncertainties and imprecision, enhancing the model’s robustness in real-world applications. This integrated approach offers a comprehensive evaluation framework, guiding informed decisions on MD-GHP project development and deployment.

• Evaluation Factor Set: The evaluation factor set is determined by selecting appropriate metrics. If there are n metrics, and the evaluation factor is u, the set is represented using Equation (39),

  $$U=\left[u_1,u_2,u_3,\cdots ,u_n\right]$$

  (39)
• Evaluation Set: For metrics with multiple rating levels, a rating scale is established based on specific criteria. Assuming there are mmm levels, the evaluation set is represented using Equation (40),

  $$V=\left[1, 2, 3,\cdots ,m\right]$$

  (40)
• Fuzzy Relationship Matrix: Given n influencing factors and mmm levels, an n × m fuzzy relationship matrix Y is constructed. Each y_i factor forms part of the matrix. The membership degree, representing the association between each metric and rating level, can be calculated using triangular membership functions using Equations (41) and (42),

  $$y_i=\left(y_{i,1},y_{i,2},y_{i,3}\cdots y_{i,m}\right)$$

  (41)

  $$Y={\left(y_1,{ y}_2, { y}_3,\cdots ,{ y}_n\right)}^T=\left[\begin{array}{ccc}{y}_{\mathrm{1,1}}& \cdots & y_{1,m}\\ ⋮& \ddots & ⋮\\ y_{n,1}& \cdots & y_{n,m}\end{array}\right]$$

  (42)
  The triangular membership function is defined by Equation (43),

  $$y_{i,k}\left(x\right)=\left\{\begin{array}{c}0, \left(if x\le a_i or x\ge c_i\right)\\ {\displaystyle \frac{x-a_i}{b_i-a_i}}, \left(if a_i<x<b_i\right)\\ {\displaystyle \frac{x-a_i}{b_i-a_i}}, \left(if b_i<x<c_i\right) \end{array}\right.$$

  (43)

  where $y_{n,m}$ with $n,m$ refers to the membership degree of the n-th influencing factor to the m-th evaluation level; $a_i$, $b_i,$ and $c_i$, respectively, represent the parameters of the triangular membership function, denoting the starting point, peak point, and ending point of the membership function.
• Matrix calculation: In fuzzy comprehensive evaluation, factors hold varying importance. By applying fuzzy operators, the membership matrix can be combined with weight results. The weighted average operator is often used, highlighting key factors while considering all metrics. The final comprehensive membership degree is calculated using Equation (44),

  $$B=W_{com}·Y$$

  (44)

  where $B$ refers to the result of the fuzzy matrix calculation; $Y$ refers to the membership degree matrix; and $W_{com}$ refers to the comprehensive matrix calculated based on EWM and AHP.
• Defuzzification and Final Evaluation Result: The final step is defuzzification, which translates the fuzzy evaluation result BBB into a crisp score or decision. This is typically accomplished using the weighted average method, which calculates a single score that represents the overall performance of the MD-GHPs. It is defined using Equation (45),

  $$S=\sum _{i=1}^m{b}_i·v_i$$

  (45)

  where $S$ is the overall score of the system; $b_i$ represents the membership degree of the system at the i-th evaluation level; and $v_i$ corresponds to the specific value assigned to the i-th evaluation level in the set $V$.

This weighted average calculation combines the influence of all factors and their corresponding membership degrees, producing a clear, quantifiable result that reflects the system’s overall performance.

4. Evaluating the Performance of Practical Projects

4.1. Analysis of the Practical Operation Performance of MD-GHPs

To validate the multistage evaluation index system for MD-GHPs developed through the AHP–EWM and fuzzy comprehensive evaluation model, the research team adopted a two-pronged approach. First, field tests were conducted on four MD-GHP projects in Shaanxi Province, China, to gather data for an EWM-based analysis, which derived objective weights. Simultaneously, a weights questionnaire was distributed to senior experts for an AHP-based analysis to obtain subjective weights. By combining the comprehensive weights with the fuzzy comprehensive evaluation model, the team analyzed and assessed the projects’ performance, effectively verifying the system’s reliability and validity.

Table 10. Results of Analysis of the Field Test Systems.

Index Layer	Unit	MG-1	MG-2	MG-3	MG-4
A1	°C	15.90	15.90	15.90	13.50
A2	°C/hm	2.76	2.85	2.85	3.01
A3	W/(m°C)	2.94	3.01	3.01	3.21
A4	kg/m3	2865.00	2990.00	2990.00	3140.00
A5	kJ/(kg·K)	910.00	860.00	860.00	940.00
A6	/	0.03	0.05	0.05	0.07
B1	kW	294.40	247.20	315.10	672.00
B2	°C/kW	0.05	0.05	0.06	0.09
B3	GJ/Year	1995.70	1946.70	2717.40	2981.10
B4	°C/Year	0.12	0.24	0.17	0.31
C1	/	4.68	5.43	6.92	6.95
C2	/	32.60	21.30	31.60	70.70
C3	/	33.80	29.30	48.80	72.70
D1	GJ/m2	0.21	0.36	0.34	0.20
D2	%	0.11	0.05	0.06	0.08
D3	%	0.18	0.12	0.09	0.13
D4	/	0.70	0.80	0.80	1.20
E1	%	0.52	0.71	0.61	0.46
E2	%	0.16	0.16	0.16	0.14
E3	%	0.00	0.00	0.00	0.00
E4	h	0.00	0.00	0.00	0.00
F1	kg/m2	19.70	30.10	20.60	11.10
F2	CNY/m2	140.00	170.00	180.00	260.00
F3	%	0.27	0.35	0.50	0.53
F4	CNY/m2	11.00	16.70	11.40	10.90
F5	Year	5.70	6.30	4.70	7.20

displays the field test results for the corresponding criterion layer indices A to F for the four MD-GHP projects. The results indicate no significant differences in resource conditions among the projects. However, the heat exchange performance of Project MG-4 was notably higher in both instantaneous and cumulative heating power compared to the others, despite its relatively poorer long-term operational performance. In terms of energy performance, Project MG-4 also showed prominent results. All projects exhibited average performance regarding building thermal demand and grid responsiveness. Although MG-4 had the lowest annual operating costs, its investment payback period was the longest. Overall, the varying performances across different dimensions highlight the necessity for further rigorous analysis based on scientific methods.

4.2. Weight Calculation of MD-GHPs

4.2.1. EWM Weight Calculation of MD-GHPs

The EWM determines objective weights by quantifying the variability of each indicator, ensuring a data-driven approach to weighting. In this study, field test results from the MD-GHP projects were used to construct a decision matrix, following Equations (20)–(22) to ensure an accurate representation of each indicator’s performance. After constructing the matrix, a standardization process was applied to normalize the data and eliminate any scale inconsistencies. Using Equations (23)–(26), the entropy for each indicator was calculated to capture the degree of variation and uncertainty within the data. These calculations culminated in the entropy weight results shown in

Table 11. Weight of Layers Based on EWM.

Index Layer	A1	A2	A3	A4	A5	A6	B1
EWM (w)	0.0415	0.0371	0.0416	0.0347	0.0550	0.0369	0.0504
Index Layer	B2	B3	B4	C1	C2	C3	D1
EWM (w)	0.0392	0.0543	0.0379	0.0463	0.0440	0.0458	0.0556
Index Layer	D2	D3	D4	E1	E2	E3	E4
EWM (w)	0.0367	0.0347	0.0359	0.0410	0.0415	0.0000	0.0000
Index Layer	F1	F2	F3	F4	F5
EWM (w)	0.0339	0.0346	0.0453	0.0399	0.0363

, which provide objective insight into the relative importance of each indicator.

Key indices emerged with high weights, including A5 at 0.055, B3 at 0.054, and D1 at 0.056. These indices emphasize the critical role of factors such as soil thermal properties, long-term heat extraction efficiency, and heating demand in evaluating system performance. The objective weighting highlights how these factors contribute significantly to the overall performance of MD-GHP systems, reflecting their importance within the operational context.

4.2.2. AHP Weight Calculation of MD-GHPs

A questionnaire was distributed to 15 experts from universities and construction and design companies to allocate subjective weights for the criterion and index layers using the AHP. The experts’ responses were used to construct judgment matrices, followed by consistency checks.

The criterion layer judgment matrix constructed based on the responses to the questionnaire is as follows:

$$A_{\mathrm{criterion}\mathrm{layer}}=\left[\begin{array}{cccccc}1& 1/5& 1/3& 3& 5& 1/6\\ 5& 1& 3& 7& 9& 1/3\\ 3& 1/3& 1& 5& 5& 1/5\\ 1/3& 1/7& 1/5& 1& 3& 1/7\\ 1/5& 1/9& 1/5& 1/3& 1& 1/9\\ 6& 3& 5& 7& 9& 1\end{array}\right]$$

The independent judgment matrices for criterion layer A to criterion layer F constructed based on responses to the questionnaire are as follows:

$$\begin{gathered} A_{\mathrm{A}\mathrm{criterion}\mathrm{layer}}=\left[\begin{array}{cccccc}1& 1/5& 1/3& 3& 5& 7\\ 5& 1& 3& 7& 8& 9\\ 3& 1/3& 1& 5& 7& 9\\ 1/3& 1/7& 1/5& 1& 2& 3\\ 1/5& 1/8& 1/7& 1/2& 1& 2\\ 1/7& 1/9& 1/9& 1/3& 1/2& 1\end{array}\right], A_{\mathrm{B}\mathrm{criterion}\mathrm{layer}}=\left[\begin{array}{cccc}1& 5& 3& 7\\ 1/5& 1& 1/3& 3\\ 1/3& 3& 1& 5\\ 1/7& 1/3& 1/5& 1\end{array}\right], \\ {A}_{\mathrm{C}\mathrm{criterion}\mathrm{layer}}=\left[\begin{array}{ccc}1& 3& 5\\ 1/3& 1& 3\\ 1/5& 1/3& 1\end{array}\right], A_{\mathrm{D}\mathrm{criterion}\mathrm{layer}}=\left[\begin{array}{cccc}1& 3& 5& 7\\ 1/3& 1& 3& 5\\ 1/5& 1/3& 1& 3\\ 1/7& 1/5& 1/3& 1\end{array}\right], \\ {A}_{\mathrm{E}\mathrm{criterion}\mathrm{layer}}=\left[\begin{array}{cccc}1& 1/5& 1/7& 1/3\\ 5& 1& 1/3& 3\\ 8& 3& 1& 5\\ 3& 1/3& 1/5& 1\end{array}\right], A_{\mathrm{F}\mathrm{criterion}\mathrm{layer}}=\left[\begin{array}{ccccc}1& 1/5& 1/7& 1/3& 1/9\\ 5& 1& 1/3& 3& 1/5\\ 7& 3& 1& 5& 1/3\\ 3& 1/3& 1/5& 1& 1/7\\ 9& 5& 3& 7& 1\end{array}\right] \end{gathered}$$

Table 12. Weight of Layers Based on AHP.

Layer	A	A1										A2										A3										A4										A5										A6
AHP (w)	0.078	0.143										0.457										0.269										0.064										0.040										0.027
Layer	B	B1															B2															B3															B4
AHP (w)	0.274	0.564															0.118															0.263															0.055
Layer	C	C1																				C2																				C3
AHP (w)	0.137	0.637																				0.258																				0.105
Layer	D	D1															D2															D3															D4
AHP (w)	0.042	0.564															0.263															0.118															0.055
Layer	E	E1															E2															E3															E4
AHP (w)	0.025	0.054															0.258															0.572															0.116
Layer	F	F1												F2												F3												F4												F5
AHP (w)	0.444	0.033												0.130												0.264												0.064												0.510

shows the calculated weights, and

Table 13. Consistency Ratio of Index Layers.

Layer	Criterion	Index
		A	B	C	D	E	F
CR	0.0771	0.0488	0.0415	0.0332	0.0415	0.0502	0.0529

illustrates the consistency ratios, all of which are below 0.1, confirming the validity of the analysis.

Based on the calculated weights, experts subjectively determined that the F Layer (Energy Saving and Economic Benefits) is the most important in the criterion layer, with a weight of 0.444. Within this category, F5 (Investment Payback Period) was identified as the most critical sub-indicator, assigned a weight of 0.510. Conversely, the E Layer indicator (Grid Dynamic Response Capacity) received a weight of only 0.025, indicating a lower level of emphasis from stakeholders regarding this sub-indicator.

4.2.3. Comprehensive Weight Calculation for MD-GHPs

Combining the results of the EWM-based analysis and the AHP-based analysis, comprehensive weights were derived (Equation (34)), integrating both subjective and objective factors. Figure 7 illustrates the comprehensive and sub-item weights for each index layer. As shown, the comprehensive weights of the B1 and F5 indices are significantly higher than those of the other indices, primarily due to the influence of subjective weighting, which elevates their perceived importance. It is evident that compared to the distribution of objective weights, the subjective weights display greater variability, likely due to differences in perception among experts from various industries.

4.3. Establishment of Evaluation Model for MD-GHPs

Building on the content of Section 2, the scoring standards for each index layer are provided in

Table 14. Scoring Criteria of Each Index Layer.

Index Layer	Unit	Poor	Normal	Good	Very Good
A1	°C	10	12.5	17.5	20
A2	°C/hm	2.12	2.41	2.85	3
A3	W/(m °C)	1.14	1.8	3.245	4.03
A4	kg/m3	1646.47	1883.95	2358.91	2596.39
A5	kJ/(kg·K)	664.13	902.28	2757.15	1616.72
A6	/	0.00265	0.0188	0.757	0.1163
B1	kW	200	300	450	500
B2	°C/kW	0.12	0.09	0.065	0.05
B3	GJ/Year	2000	2500	3500	4000
B4	°C/Year	0.8	0.6	0.3	0.2
C1	/	5	5.5	6.5	7
C2	/	40	50	70	80
C3	/	40	50	70	80
D1	GJ/m2	0.3	0.25	0.16	0.12
D2	%	20%	17.50%	10%	5%
D3	%	5%	4.50%	3.50%	3%
D4	/	1.5	1.25	0.75	0.5
E1	%	40%	50%	70%	80%
E2	%	10%	15%	25%	30%
E3	%	20%	25%	35%	40%
E4	h	2	3	5	6
F1	kg/m2	30	25	15	10
F2	CNY/m2	300	250	150	100
F3	%	20%	25%	35%	40%
F4	CNY/m2	15	12.5	7.5	5
F5	Year	8	6.5	4	3

. For each indicator, after determining specific evaluation levels, reference is made to set values from existing standards or to quartiles derived from nationwide data distribution. Depending on the characteristics of each indicator, appropriate mathematical and statistical methods, such as equal spacing interpolation or equal proportion floating, are used to determine the remaining levels. These scoring standards form the basis for determining membership degrees using the triangular fuzzy function method.

Based on Equations (32) and (33), during the fuzzy composition process of the factor sets, the hierarchical weight of each index is determined using the AHP–EWM approach. The comprehensive membership degree of the factor set in the criterion layer represents the overall fuzzy evaluation result. After processing the data from the decision set, a comprehensive score is obtained, which reflects the actual performance of the system. This comprehensive score is a useful tool for identifying weak points in the operation of MD-GHP systems. Figure 8 illustrates the evaluation results of the multistage index system for the four field test projects. From the figure, it is clear that the C and E layers are relatively weak, with most of the sub-indices falling in the mid-level range, and no sub-indicator achieving exceptional results.

4.4. Analysis of Evaluation Results for MD-GHPs

Based on the model and evaluation results, this section provides an in-depth analysis of the factors affecting performance scores across the four MD-GHP projects. Figure 9 illustrates the comprehensive scores for each project: MG-1 scored 61.56, MG-2 scored 58.33, MG-3 scored 72.73, and MG-4 scored 78.41. In this figure, green represents sub-indices rated as “Very Good,” while red denotes “Poor” ratings. Across the evaluation categories, several key factors contribute to performance differences.

In terms of resource conditions (A-Layer), all projects performed well due to thorough planning and site selection, ensuring stability for effective system operation. For heat transfer performance of the MDBHEs (B-Layer), MG-3 and MG-4 demonstrate superior scores, particularly in index layer B3, attributed to factors such as lower inlet water temperature, higher flow rates, and greater DBHE depths. MG-4, in particular, benefits from its DBHEs reaching 2800 m in depth and an intermittent operational mode (11 h on, 13 h off), which promotes greater heat extraction and higher outlet water temperature. Through index layer B4, it can also be seen that all four projects have the potential for long-term sustainable operation. The energy performance of the MD-GHPs (C-Layer) also highlights MG-4 as an outlier, as it employs a specially designed, second-generation heat pump capable of handling large temperature variations (20 to 35 °C) and partial load ratios, unlike the conventional SD-GHPs used in MG-1 to MG-3. Additionally, MG-4’s transmission and distribution system benefits from an optimized design, enhancing temperature differentials and reducing water resistance. Regarding building space heating demand (D-Layer) and grid dynamic response capacity (E-Layer), all projects reveal improvement opportunities related to building design, heating pipe network enhancements, and operational strategies to align with grid demands. Lastly, energy-saving and economic benefits (F-Layer) show MG-4 outperforming the others due to its use as an office building, where separate heating systems in traditional offices generally lead to higher energy consumption. The use of MD-GHPs in office settings, therefore, demonstrates notable economic and energy-saving advantages.

In summary, these findings underscore the major factors impacting disparities in project performance, including system depth, operational modes, heat pump design, and tailored operational strategies. Future optimizations should focus on these elements to improve MD-GHP system performance across multiple metrics.

4.5. Analysis of Energy-Saving and Emission Reduction of Different Heating Sources

To better evaluate the energy saving and CO₂ emission reduction effect of MD-GHPs, the energy performance of conventional heat pumps and optimized heat pumps in MD-GHPs have been compared with heat pumps in ASHPs, SD-GHPs, and gas boilers.

Table 15. Analysis of energy cost, primary energy consumption, and CO₂ emissions.

	MD-GHPs	ASHPs	SD-GHPs	Gas Boilers
Designed heating capacity (kW)	7560
Cumulative Qc (GWh)	17.48
Energy efficiency of System	5.09	2.40 [28]	2.81 [67]	50 for user-side water pumps [28]
Electric consumption (GWh)	3.43	7.28	6.22	0.35
Gas consumption (Million Nm3)	/	/	/	1.84
Energy Cost (Million CNY)	2.75	5.83	4.98	6.90
Primary Energy (tons of standard coal equivalent)	1064.60	2257.83	1928.40	2555.58
CO2 Emissions (tons)	2060.51	4370.00	3732.38	4000.16
Initial costs (Million CNY)	26.05	7.56 [28]	13.61 [68]	6.35 [28]
Static incremental payback period (Year)	/	6.0	5.6	4.7

shows a comparison of the energy cost, primary energy consumption, and CO₂ emissions of these systems. For the analysis, the cumulative heating consumption and the energy efficiency of MD-GHPs in MG3 have been applied, which reach 17.48 GWh and 5.09, respectively.

For the initial cost of MD-GHPs, based on feedback from the project owner, the initial cost per MDBHE reached about CNY 2.0 million per MDBHE with a depth of 2500 m in Xi’an. Then, considering the device and installation costs of MDBHEs and heat pump systems, the total initial cost of MD-GHPs in MG3 reached CNY 26.05 million. It can be seen that the energy efficiency of MD-GHPs was superior to that of the other three types of heat pumps. Thus, the electric consumption reached only 3.43 GWh; while the electric consumption of ASHPs, SD-GHPs, and user-side water pumps in gas boiler systems reached 7.28 GWh, 6.22 GWh, and 0.35 GWh, respectively. Notably, the electric consumption of MD-GHPs is 52.8% and 44.8% lower than ASHPs and SD-GHPs, respectively. Simultaneously, gas boilers require a gas consumption of 1.84 million Nm³. In terms of primary energy consumption, the MD-GHPs produced energy saving rates that were 52.8%, 44.8%, and 58.3% higher than those of ASHPs, SD-GHPs, and gas boilers, respectively. In addition, the CO₂ emission reduction rate of MD-GHPs was 52.8%, 44.8%, and 48.5% higher than those of ASHPs, SD-GHPs, and gas boilers, respectively.

From the perspective of operational energy cost, taking an electricity price of 0.8 CNY/kWh and a gas price of 3.6 CNY/Nm³ in Xian, the operational energy cost of MD-GHPs reached CNY 2.75 million while the energy cost of ASHPs, SD-GHPs, and gas boilers reached CNY 5.83 million, CNY 4.98 million, and CNY 6.90 million, respectively. Thus, the energy cost of MD-GHPs was 52.8%, 44.8%, and 60.2% lower, respectively. To further evaluate the economic effect, the initial costs of different heat sources were assumed according to a previous study [28] at CNY 26.05, 7.56, 13.61, and 6.35 million for MD-GHPs, ASHPs, SD-GHPs, and gas boilers, respectively. Although the initial costs of MD-GHPs were significantly higher, due to the higher energy efficiency and significant energy saving effects, the static incremental payback period of MD-GHPs was about 6.0, 5.6, and 4.7 years compared with ASHPs, SD-GHPs, and gas boilers, respectively. That is to say, MD-GHPs are superior to other heat sources in terms of energy cost, primary energy consumption, and CO₂ emissions, making them more suitable for buildings space heating from the perspective of long-term operation and development.

5. Conclusions

This study developed a comprehensive, multistage evaluation index system tailored for MD-GHPs, addressing a crucial gap in current research. Given the significant energy consumption and environmental impact associated with space heating in China, MD-GHPs present a viable solution for achieving low-carbon and sustainable heating. The proposed evaluation index system offers a refined approach to identifying weaknesses in MD-GHP applications, facilitating a more scientific and practical analysis for system design, implementation, and operation. The key findings are as follows:

(1)

Based on the specific characteristics of MD-GHPs, a multi-level evaluation index system was constructed by referencing relevant standards and literature, and by extensive field tests conducted by the research team. This system encompasses six core dimensions: geothermal resources, the heat transfer performance of the MDBHEs, the energy efficiency of the heat pump systems, building space heating demand, grid dynamic response capability, and energy-saving and economic benefits. This robust framework comprises 26 specific indices, effectively evaluating both the heat-source and user sides of the system under real-world conditions. Additional assessments include the quantity of clean energy utilized, energy resilience, and economic benefit evaluations, enhancing the quantitative accuracy of the system’s performance assessment.

(2)

Using both expert survey results and case study data, subjective and objective weighting methods were applied to determine indicator weight coefficients. A linear recombination method was employed to harmonize the strengths of both weighting approaches, mitigating the limitations of each. The results show that among the six dimensions, “energy-saving and economic benefits” has the highest weight (0.317). The top three indices within the index layer are the system investment payback period (0.131), system energy-saving rate (0.081), and comprehensive heat transfer resistance (0.036).

(3)

A fuzzy comprehensive evaluation model was established to integrate these weighted indices into a multi-criteria evaluation of MD-GHP performance. This model enables an overall assessment of system performance during building operations, highlights technical weaknesses, and offers insights for optimizing system operations and guiding investment decisions.

(4)

Evaluation of four MD-GHP projects based on real-world data revealed the following performance scores: MG-1 scored 61.56 (Good), MG-2 scored 58.33 (Normal), MG-3 scored 72.73 (Good), and MG-4 scored 78.41 (Good). These evaluations indicated that the heating system performance (C-level) and grid dynamic response capability (E-level) were the weakest areas, necessitating improvement. Future MD-GHP projects should enhance heating system performance through retro-commissioning (RCx), operational strategies, and intelligent control mechanisms. Furthermore, optimizing power grid dynamic response and system resilience will be crucial for developing more efficient, green, and low-carbon MD-GHP systems.