A numerical study on performance efficiency of a low-temperature horizontal ground-source heat pump system

Exploitation of shallow ground and its low-grade heat potential is fundamental to designing 5th generation district heating and cooling (DHC) networks. Horizontal ground-source heat pump (HGSHP) systems are a common way to utilize shallow geothermal energy. Realistic estimation and prediction of performance of a HGSHP system and shallow ground thermal behaviour should consider the whole system including building heating and cooling load, heat pump and ground heat exchanger, and the ground. This should be accompanied by realistic atmospheric and ground conditions. In this paper, a three-dimensional coupled thermal – hydraulic model with realistic boundary conditions adopting a whole system approach is presented. Dynamic heat pump coefficient of performance (COP) that depends on seasonal variation of heating/cooling demand and ground conditions are also considered. Model validations are conducted against experimental and analytical results in literatures. The model is applied for evaluating a HGSHP system to support development of a 5th generation DHC network on a potential site in the UK. Several influencing factors, such as ground moisture transfer, building thermal load mode, buried depth of ground loops, and initial ground temperature profile are studied to assess performance efficiency of the HGSHP system and evolution of ground thermal behaviour in response to heat extraction or rejection into the ground. The results show that 5% of the monthly total heat demand of the site could be met by the designed HGSHP system, consisting of 200 U-shaped ground loops buried at the depth of 3 m and pure water as the heat carrier. Overlooking the ground moisture transfer or hyperbolizing the ground saturation would overestimate the load-carrying capacity of the HGSHP system. The HGSHP system is more efficient with a higher heat pump COP under the heating and cooling mode than under the heating-only mode. Predicted performance of the HGSHP system improves with buried depth of the ground loops. The results also show that a 1 ℃ increase in the undisturbed ground temperature could suffice up to 8% of the monthly total heat demand of the site.


Introduction
Ground-source heat pump (GSHP) systems that utilise shallow ground as a heat source or sink offer energy-efficient heating and cooling solutions as well as minimise environmental impact by reducing greenhouse gas emission [49,12,45,39].Exploitation of shallow ground and its low-grade heat potential is fundamental to designing 5th generation district heating and cooling (DHC) networks.A closed-loop, horizontal GSHP or HGSHP system exerts a good balance between efficiency and cost [45,39].Thus, a comprehensive performance analysis of such system is conducted within the scope of this study.
Performance efficiency of the HGSHP system and thermal behaviour of the shallow ground depend on a number of factors, which can be broadly classified into four categories: i) site conditions, such as the local climatic conditions, ground temperature variation, and soil saturation/ moisture change [32,51,37,30]; ii) ground characteristics, such as soil thermal conductivity, density, and specific heat capacity, which vary greatly between different soils [41,10,56,47], iii) system configurations and design, such as pipe length, diameter, spacing, buried depth, the thermal load of the building, and integration with other renewable energy sources, such as solar energy, which have considerable effects on the performance of the HGSHP system [6,[17][18][19]27,28,36,52]; and iv) heat pump coefficient of performance or COP, which regulates the thermal load on ground loop systems [55,45,52].
Experiments have been conducted to study influential factors on the performance efficiency of the HGSHP system.For example [17,18], constructed a ground-source heat pump system coupled with two horizontal ground heat exchanger circuits, which were buried at different depths below the ground surface.The experiment results showed the system's COP and the heat pump's COP of the system coupled with the 2m-deep horizontal ground heat exchanger circuit were 2.82 and 3.42, respectively, both higher than those (2.68 and 3.13) of the system coupled with the 1-m-deep horizontal ground heat exchanger circuit.In addition, the energy efficiency and the exergetic efficiency of the system coupled with the 2-m-deep horizontal ground heat exchanger circuit were 2.8 and 56.3%, respectively, also both greater than those (2.5 and 53.1%) of the system coupled with the 1-m-deep horizontal ground heat exchanger circuit [19].Furthermore, using the measured experimental data, such as air temperatures, fluid temperatures, and ground temperatures as input, the performance of the horizontal ground source heat pump was investigated and predicted through artificial neural network approaches [21][22][23][24][25][26].
Numerical investigations on the efficiency of HGSHP systems in relation to heat extraction and injection of geothermal energy are reported in literatures.Esen et al. [20] developed a two-dimensional (2D) ground heat transfer model for determining the temperature distribution near horizontal pipes.They observed that the temperature distribution surrounding pipes influences the heating load of the HGSHP system and its COP.However, influences of the heat transfer in the direction parallel to the pipes, rainfall, and true atmospheric boundary conditions on the system behaviour were ignored.Pulat et al. [49] also built a 2D numerical model to study the thermal interaction among the horizontal pipes.However, the surface temperatures of pipes between inlet and outlet pipes were fixed owing to the 2D steady-state assumption and their values were estimated based upon measured inlet and outlet water temperatures, which contradicted the dynamic fluid temperature as fluid circulated and exchanged heat with the adjacent ground.Therefore, compared to 2D modelling, three-dimensional (3D) modelling is more suitable and realistic for studying the performance of a HGSHP system and the thermal behaviour of shallow ground.
3D simulations have been conducted by researchers to study the influential factors on the thermal performance of the HGSHP system and to validate design methods.Dasare and Saha [14] developed a 3D numerical model to predict the thermal performance of three configurations of horizontal ground loops (linear, helical and slinky) under the influences of buried depth, fluid flow rate, and soil thermal conductivity.They found the soil thermal conductivity is the most vital parameter in the heat transfer process and the helical geometry performs better than the other two geometries.A 3D numerical model was established by Li et al. [45] to evaluate the operation characteristics of the horizontal spiral-coil GSHP system by the effects of subsurface factors, daily variations of load, and operation models.They concluded that soil thermal conductivity and pipe spacing are two most important subsurface factors, daily variations of load can be neglected, and continuous operation provides the best performance.Kim et al. [39] proposed a design equation for the design length of the horizontal spiral-coil GSHP system, and a 3D computational fluid dynamic (CFD) simulation was performed to validate the applicability of the proposed equation.For the above thermal modelling and design of the GSHP system [14,45,39], the analytical expression, which is a sinusoidal equation, was adopted to approximately represent the ground temperature profile, namely, the variation of ground temperature with depth and time.The expression was analytically derived by Kusuda and Achenbach [40] for a semiinfinite solid considering heat conduction below the ground surface because of the ambient temperature.However, the ground temperature profile is highly dependent on the local climatic conditions [51,37].The other climatic conditions, such as wind speed, rainfall, and short-wave radiation cannot be reflected by the analytical expression.Furthermore, parameters in the analytical equation are not easy to determine, while the climatic variables are commonly monitored in the field and representative values can be obtained for most regions globally from metrological data.
Existing models on HGSHP systems often simplify ground characteristics.For example, Esen et al. [20] ignored soil moisture transfer processes to obtain the temperature distribution in the vicinity of pipes.Pulat et al. [49] assumed constant soil conductivity in their numerical model when analyzing thermal interactions among the horizontal pipes.The design length of the horizontal spiral-coil GSHP system proposed by Kim et al. [39] was derived from ground heat transfer alone.In addition, Kayaci and Demir [38] and Sedaghat et al. [52] adopted constant thermal properties of soils and only considered the heat transfer process in modelling long-term performance analysis of the HGSHP system.However, soil is a three-phase system consisting of solid aggregates, pore-water, and pore-air.In reality, heat and moisture exchanges occur constantly through the soil-atmosphere interface and soil thermal behaviour depends on responses of these constituents, both individually and collectively [30].Moreover, experiments have revealed that the thermal properties of soil are influenced by its water content.For example, the thermal conductivity of sandy soil can be reduced from 2.65 W/m/K in saturated condition to 0.90 W/m/K in dry soil [2].The thermal conductivity of Basaltic Clay (silty clay) was decreased from 1.34 W/m/K to 0.32 W/m/K when the saturation dropped from 100% to 13% [3].Therefore, moisture transfer and dynamic soil thermal properties should be considered for true ground representation and accurate evaluation of long-term performance of a HGSHP system as well as thermal potentials of shallow ground, which is not sufficient in literatures.
Along with focusing on the features of horizontal ground loops, such as configurations or the working characteristics of a HGSHP system, the heat pump COP significantly affects the performance efficiency of the HGSHP system and the thermal behaviour of the shallow ground.For instance, Li et al. [45] compared the inlet and outlet fluid temperatures in the horizontal ground loops with and without the heat pump.They concluded that the inlet and outlet fluid temperatures with the heat pump are higher in the cooling mode but lower in the heating mode.Since the heat pump is an important part of the HGSHP system, numerical models should consider the heat pump COP and its operational variations when analysing the long-term performance of the whole system.
In this paper, a 3D coupled thermal-hydraulic (TH) model with realistic boundary conditions is presented.The aim is to investigate performance efficiency of a HGSHP system and thermal response of a shallow ground under both heating and cooling modes.The model takes into account the whole system, broadly containing the modules, such as heat and moisture transfer into and within a shallow ground (module 1), heat transfer within the horizontal ground loops (module 2), and the coupling of the horizontal ground loops to a heat pump (module 3).The proposed model is significantly advanced than the existing models since it explicitly represents ground thermal and moisture transport processes and their combined influence on the HGSHP system performance.In addition, this model considers the change of thermal properties of soil with its types and variations in the ratios of solids, pore water and pore air.Accurate design and assessments of low temperature, groundsourced, 5th generation DHC networks closely depend on detailed ground energy, mass transfer processes and coupled atmospheric conditions.However, numerical models those include the aforementioned features and processes are rare in literatures.
In the application of the proposed model variation of soil layers along the depth of a shallow ground, selected for potential development of a 5th generation DHC network in the Warwickshire County, UK, is taken into consideration.Local atmospheric data, such as, solar radiation, rainfall, humidity, air temperature, and wind velocity are taken into account in the model simulations.Along with the soil-atmosphere boundary, processes at the soil-soil boundary such as heat and moisture transfer from the areas adjacent to the site is considered for realistic representation of in-situ conditions.Moreover, the coupling of the ground loops and the heat pump is achieved by the heat pump COP model proposed by Staffell et al. [55].

W. Gao et al.
The structure of this paper is as follows.The theoretical and numerical formulations related to moisture and heat transfer within unsaturated soils, ground surface boundary, ground-loop boundary, and heat pump COP were presented first.Then the numerical solutions were explained in detail.Following that the model was validated against experimental and analytical results from literatures.Finally, the proposed model was used to investigate the effects of various influential factors, including ground moisture transfer, building thermal load mode, buried depth of ground loops, and initial undisturbed ground temperature profile, on the performance efficiency of a HGSHP system and the thermal behaviour of the shallow ground due to the extraction and injection of geothermal energy.

Theoretical and numerical formulations
Fig. 1 illustrates the horizontal ground-source heat pump (HGSHP) system considered in this study.It consists of a shallow ground, horizontal ground loops parallelly buried in the ground, a heat pump, and a heat distribution unit in a building.The general assumptions made in this paper to develop the numerical model are as follows: soil is homogeneous and isotropic, ground loop pipes are equally spaced, and the distance between pipes is big enough to avoid thermal interference with each other.

Moisture and heat transfer within unsaturated soils
Moisture flow within unsaturated soils is described as a two-phase process, compromising of both liquid and vapour flows.The general equation for the moisture flow can be described as: where θ l is the volumetric water content, θ a is the volumetric air content, t is the time, ∇ is the gradient operator, ρ l is the density of the water, and ρ v is the density of vapour.The velocity of water v l is calculated by Darcy's law: in which u l is the pore-water pressure, γ l is the unit weight of water, y is the elevation, and K l is the unsaturated hydraulic conductivity, which is modelled using the Brooks and Corey [8] Model herein: where K ls is the saturated hydraulic conductivity, θ ls is the saturated water content, and η is the shape parameter.The Van Genuchten [59] Model is used to characterize the Soil Water Characteristic Curve (SWCC) of soils, and there is: where h is the pressure head, h d is the scale parameter, and n is the shape parameter.
The velocity of vapour v v is described by the following equation [48]: in which D atms is the molecular diffusivity of vapour through air, v v is a mass flow factor, τ v is a tortuosity factor, and ∇ρ v is the spatial vapour density gradient.Heat transfer within unsaturated soils is considered mainly via conduction and convection.The temporal derivative of the heat content is equal to the spatial derivative of the heat flux, and there is: On the left hand of Equation ( 6), L is the latent heat of vaporisation, and ϕ is the soil porosity.H c is the heat capacity of unsaturated soil at a reference temperature T r , which is calculated as follows [29]: where C ps , C pl , and C pv are the specific heat capacities of the solid, water, and vapour, respectively, S l and S a are the degrees of saturation of water and pore air, respectively, and ρ s is the density of the solid.
On the right hand of Equation ( 6), the first term describes the thermal conduction in accordance with Fourier's law, the second term describes the latent heat flow associated with vapour movement, and the third term describes the heat convection in terms of the liquid phase movement and vapour phase associated with a vapour pressure gradient [57].λ T is the thermal conductivity for the unsaturated soil, which can be obtained by the soil's components as the following model [58]: Fig. 1.Horizontal ground-source heat pump system considered in this study and essential parts of the system.
where λ s , λ w , and λ a is the thermal conductivity corresponding to the solid, water, and air, respectively, and χ s , χ w , and χ a is the volume fraction corresponding to the solid, water, and air, respectively, which can be presented by the following expressions:

Boundary conditions
Fig. 2 shows the representative cross-section in the ground of the HGSHP system in Fig. 1.The horizontal parallel pipes are placed with a spacing S in the ground at the same burial depth H.The length, outer radius, and thickness of pipes are L, R, and b respectively.The fluid flow rate r f is assumed to be the same in each pipe.
As can be seen from Fig. 2, there are two boundary conditions required for the modelling of the ground behaviour in response to geothermal energy extraction or injection via the HGSHP system: i) the ground surface boundary, and ii) the ground-loop boundary.

Ground surface boundary
The ground surface boundary involving heat and moisture exchange between the ground and atmosphere needs to be considered due to the large area covered and the shallow buried depth of horizontal ground loops.
At the ground surface, the energy balance equation can be developed as follows [9,62]: where H NE is the net radiant energy flux absorbed or emitted at the ground surface, H Absorbed SW is the absorbed shortwave radiation flux at the ground surface, H Absorbed LW is the longwave radiation flux absorbed at the ground surface, H Emitted LW is the longwave radiation flux emitted from the ground surface, H SEN is the sensible heat radiation flux at the ground surface, and H LE the latent heat radiation flux at the ground surface.
The equation of H Absorbed SW can be expressed as below [63]: where ε SW is the shortwave reflection factor associated with the ground surface, and it is taken as 0.215 herein, and H SW is the shortwave solar radiation.
The value of H Absorbed LW can be obtained by the following equation [64,65]: in which σ is the Stefan-Boltzmann constant, C cloud is the cloud cover coefficient, and it is assumed that 1 +0.17C 2 cloud = 1.1 herein, T a is the ambient air temperature adjacent to the ground surface, and ε a LW is the long-wave emissivity of the air at ground level, which equals to 9.2 × 10 − 6 • T 2 air .Based on the Stefan-Boltzmann's law [66], the value of H Emitted LW is calculated below: where ε LW is the long-wave emissivity of the ground, and it is taken as 0.96 herein.
For H SEN and H LE , the following relations are used, respectively [63]: where ρ a is the air density, C pa is the specific heat capacity of air, u y is the wind speed, and C y is a drag coefficient and taken as 0.016 herein, L is the latent heat of vaporization, and E is the evaporation flux (saturated state), which can be calculated as follows [63,67]: in which q is the specific humidity of the soil at the ground surface, and q a is the specific humidity of air.As the ground surface would enter the unsaturated state, the following modification is made to Equation ( 16) and there is [68,69]: where E AE is the actual evaporation flux (unsaturated state), h g is the relative humidity of the ground surface, and h a the relative humidity of air at the ground surface.
Assuming the moisture flux at the ground surface is a hydrological process, the moisture balance at the ground surface can be presented below [70]: where M NM is the net moisture flux at the ground surface, P is the rainfall, and R RO is the run-off.

Ground-loop boundary
The theoretical formulation of the ground-loop boundary contains two parts: firstly, the fluid temperature profile along the ground-loop Fig. 2. Schematic showing boundary conditions in the representative cross-section in the ground of the HGSHP system in Fig. 1.
W. Gao et al. length is predicted based on the conservation of energy; secondly, the heat flux between the fluid and the adjacent ground is calculated.
As shown in Fig. 3, taking a U-shaped ground-loop as an example, its two legs are discretised into a series of uniform control volumes of length dL.The connecting part of the two legs is not considered due to its small size.The fluid temperature is assumed to be uniform across each control volume.Due to small pipe diameter, high convective heat transfer coefficient, and turbulent flow inside the pipe, the fluid temperature is assumed to be uniform in a section perpendicular to the pipe axis [32,5,45,42,38].Therefore, the temperatures of the control volumes constitute the fluid temperature profile, which changes along the ground-loop length.
For the first and last control volumes in the ground-loop, the fluid temperature corresponds to the inlet fluid temperature T f,i and the outlet fluid temperature T f ,o , respectively.Within a time period dt, geothermal energy would be removed from the circulating fluid through the heat pump in heating mode or added to the circulating fluid in cooling mode.The following relationship between inlet fluid temperature and outlet fluid temperature can be obtained: where Q GL is the thermal load of the ground-loop (positive in heating mode, negative in cooling mode), C pf is the specific heat capacity of the fluid, and ρ f is the fluid density.
For the jth control volume in the ground-loop, the fluid temperature change for a given time period dt can be calculated as below: in which Q j is the heat flow rate between the jth control volume and the adjacent ground, which can be expressed: where R res is the thermal resistance of pipes, and λ p is the thermal conductivity of pipes.Based on Fourier's law, the heat flux at the ground-loop boundary can be calculated by the following equation: The application of Equations ( 19) and ( 20) allows the updated fluid temperature in each control volume to be calculated.Based on the updated fluid temperature profile and Equation ( 22), the heat flux will then be calculated and used to model the ground thermal behaviour.

COP of heat pump
By fitting the industry-average data of the heat pump COP and operation temperatures for UK conditions, Staffell et al. [55] proposed an empirical model for predicting the heat pump COP in heating and cooling modes, respectively, which is adopted in the literature [31,52] and this study.The heat pump COP can be obtained by the following equations: where ΔT is the temperature difference, T hot is the temperature of the supplied hot water to the building, and T chilled is the temperature of the supplied chilled water to the building.Based on the literature [50],T hot usually varies from 45 • C to 40 • C in winter and T chilled usually ranges from 7 • C to 12 • C in summer.In this study, the values of T hot and T chilled are taken as 42.5 • C and 9.5 • C, respectively.
Given the building thermal load Q Building and the heat pump COP under heating or cooling mode, the thermal load of the ground loop Q GL can be obtained as follows [50,52]: Through Equation ( 23), the ground loops and the heat pump are coupled via the temperature difference between the cold and hot sources.From Equations (24a) and (24b), the thermal load of the ground loop is lower than that of the building in heating mode, and vice versa in cooling mode.

Numerical solutions
The coupled TH model presented in this paper has been developed within the framework of a bespoke thermo-hydraulic-chemical-mechanical (THCM) modelling code, namely, COMPASS (COde of Modelling PArtially Saturated Soils) [57,30].In the model, the governing equations are expressed in terms of the primary variables, i.e., pore-water pressure u l and temperature T. Equation ( 1) and Equation ( 6) can be expressed as follows: Fig. 3. Schematic showing the discretisation of the ground-loop and fluid temperature profile (not scaled).
W. Gao et al.

C TT ∂T ∂t
where C and K terms represent storage and flux, respectively, and detailed in the Gao et al. [30].
Fig. 4 illustrates the schematic flowchart of the numerical solution.As shown in the figure, a multiple time-steps algorithm, including the global time step loop and the local time step loop, is adopted to obtain the thermal behaviour of the ground and the fluid temperature profile.Within the framework of COMPASS, the Galerkin finite-element method [61] is adopted to spatially discretize the governing equations of the coupled TH model, namely, Equations ( 25) and (26), and an implicit mid-interval backward difference time-stepping algorithm is employed for temporal discretisation.The discretised system of linear equations is solved iteratively using a predictor-corrector algorithm [15] to obtain the ground thermal behaviour at the global time-step level, including the ground temperature and pore-water pressure distributions.Between two consecutive global time-steps, the local time-step is prescribed, whose length equals the global time-step, the fluid temperature profile is calculated under certain fluid flow rate, pipe configurations and properties, fluid properties, and building thermal load.By repeating this multiple time-steps procedure, the fluid circulation within the groundloop can be mimic.Meanwhile, the updated fluid temperature profile is used to update the heat flux at the ground-loop boundary to obtain the ground thermal behaviour.

Model validation
In the previous research [30], the validation of the coupled TH model with the ground surface boundary was carried out by comparing the simulated evaporation at the soil atmospheric interface against the measurements of an in-situ experimental study conducted in Southern France by Enrique et al. [16] and Calvet et al. [11].Two more cases in literatures [35,44] were employed to validate the proposed model with the ground-loop boundary herein.
4.1.On-site heating experiment with a vertical borehole [35] Hu [35] conducted a heating experiment using a U-shape pipe.Information on geologic conditions, soil physical parameters, and borehole and fluid parameters were well documented and listed in Table 1 and Table 2, respectively.A vertical borehole was drilled in the ground consists of three layers of soils, e.g., backfill soil, clay, and fine sand with the thickness of 20 m, 18 m, and 25 m, respectively.Saturated soil condition prevailed due to the presence of groundwater table at the ground surface, and no groundwater flow was monitored during the experiment [35].A U-shape pipe was installed inside the 0.07-m-radius borehole that filled by grout material.
Due to the lack of information on the physical parameters of grout material, the average values of the three layers of soils were used in the study.The grout material's effective thermal conductivity was taken as Fig. 4. Flowchart of the numerical solution.
W. Gao et al.
1.50 W/m/K and the initial temperature of the soils was set as 15.5 ℃ (288.65 K), which were as identical as in the simulation by Hu [35].In addition, heat transfer in the grout material was regarded as being purely conductive.Fig. 5 compares the outlet fluid temperature predicted by the presented model with the experimental and analytical results of Hu [35].As can be seen from the figure, the simulated results obtained by the present model were in good agreement with Hu's modelling results, and the difference was less than 0.02 ℃.In contrast to the experimental data, both the present model and Hu's modelling showed deviations, which perhaps were associated with the grout material properties and heat transfer mechanism.The difference decreased with time and was only 0.2 ℃ at the end of the experiment.

4.2.
Laboratory heating experiment with a horizontal ground heat exchanger [44] Li et al. [44] built an indoor experimental apparatus.The rectangular cuboid apparatus was 6.25 m × 1.5 m × 1.0 m (length × height × width) and installed laterally.A copper U-shape pipe was installed and penetrated two soil layers horizontally, e.g., 3-m thick soil and 3.25-m-thick clay.Information on soil physical parameters, and horizontal ground heat exchanger and fluid parameters were given in Table 3 and Table 4, respectively.The initial soil temperature was 12.23 ℃.In the experiment, hot water was circulated at a constant flow rate (0.64 m/s) in the horizontal ground heat exchanger, and the inlet and outlet fluid temperatures were monitored for 48 h.
The experimental apparatus was enveloped with 240 mm-thick brick walls and a 30 mm-thick thermal insulation mortar was added to the

Table 2
Basic parameters of the borehole and U-shape pipe [35].

Table 3
Physical parameters of soils [44].interior [43].Therefore, it is assumed in the study that the surface temperature of the apparatus remained constant at 12.23 ℃.
The outlet fluid temperature predicted by the presented model and the experimental results obtained by Li et al. [44] are compared in Fig. 6.As shown in the figure, the simulated results from the proposed model and the experimental data from Li et al. [44] agreed well.The maximum temperature difference was less than 0.35 ℃, and it decreased with time.
Based on the model validations, the good agreement between the predicted and observed results indicate successful implementation of the proposed model in COMPASS code, and the accuracy of the proposed model in predicting heat transfer processes of ground-source heat pump system in complex geological conditions.

Applications
The proposed model was used to study the performance efficiency of a HGSHP system and the thermal behaviour of the ground due to the heating and cooling process.The application testing site is located within the University of Warwick campus in the Warwickshire County, UK.The site is selected for a potential development of a 5th generation DHC network under the "Integrated Heating and Cooling Networks with Heat-sharing-enabled Smart Prosumers" project (Grant number: EP/ T022795/1).The ground available for the development is around 100000 m 2 for this site.An initial estimation suggested installing horizontal U-shaped ground loops with a total number of 200 would be ideal for the site.The heating demand data was obtained from the University's estate department.The local atmospheric data was obtained from the meteorological station at Church Lawford [46].The ground profile data were obtained from the British Geological Survey [7] borehole database.Three borehole logs at the vicinity of the site categorized the ground into three layers, e.g., Layer 1: 0-0.3 m sandy clay loam, Layer 2: 0.3-2.4m silty clay, and Layer 3: greater than 2.4 m mudstone.
The system consisted of 200 U-shaped ground loops that were buried in parallel at the same depth (3.0 or 3.5 m).The length and spacing of the pipes of each U-shaped ground-loop were 200 m and 1 m, respectively.The pipe outer radius was taken as 0.02 m, the most common pipe size used in the UK [13].In addition to that, pure water was chosen as the fluid in pipes due to its non-toxicity, low expense, and good thermal properties [4,53].Based on the cold running water in the UK and the ground temperature, the initial fluid temperature was taken as 12 ℃.The fluid flow rate (2.5 m 3 /h) was determined based upon a site experiment conducted by Hepburn et al. [33].

Model domain
Due to the periodic configuration of the horizontal ground heat exchangers, one representative U-shaped ground-loop was simulated to reduce the computational cost.As shown in Fig. 7, simulations were conducted on a domain with the size of 2.0 m × 4.0 m × 200.0 m (2S × D × L).From the figure, the Cartesian coordinate system was placed at the ground surface and the middle of the domain width.Considering the balance between the computation efficiency and the accuracy of results, Fig. 7 shows the optimal mesh after conducting simulations with different mesh sizes.The discretised domain constitutes 31,300 hexahedral elements which are connected by 35,190 nodes.In the x-y plane, the domain area around the pipes was discretised into small elements (0.01 ~ 0.02 m) to consider the active heat exchange between pipes and the adjacent ground.Along the z direction, the domain length was divided into 50 segments with a uniform length of 4 m, which was also

Table 4
Basic parameters of the copper U-shape pipe and fluid [43,44].the length of control volumes of the pipes.In the simulations, fluid is assumed to enter from the right pipe and exit from the left pipe.As shown in Fig. 7, the U-shaped ground loop was 3.0 m below the ground surface, and an observation point P was selected at the location (0.52, − 3.0, 0) to monitor the ground temperature.The observation point P is right next to the right pipe where inlet flow enters, therefore, the ground temperature at the observation point P can represent the lowest ground temperature during the heating process.For the case with a buried depth of 3.5 m, the corresponding location was (0.52, − 3.5, 0).Table 5 lists the physical parameters of the three soil layers at the application testing site, where the parameter values were determined from data in the literature [41,10,54,47].The parameters of ground-loop pipes and fluid are listed in Table 6 and Table 7, respectively.Based on the fluid properties, flow rate and pipe size in Tables 6 and 7, a turbulent flow would be established in the horizontal pipe with a Reynolds number of approximately 12,400 and subsequently a sufficient convective heat transfer between the fluid and the pipes.Additionally, the values of thermal conductivity for pore-water and pore air were taken as 0.57 and 0.025 W/m/K, respectively.The value of specific heat capacity of vapour was taken to be 1870 J/kg/K.

Initial and boundary conditions
As shown in Fig. 7, the symmetry boundary condition was applied to the left and right sides of the domain.On the surface of the domain,  atmospheric data in Fig. 8 [46] were prescribed.The mean undisturbed ground temperature at 4 m below the ground surface was about 12 ℃ [71], therefore, a fixed temperature boundary was applied at the bottom of the domain to consider the undisturbed ground temperature.In addition, the pore-water pressure at the bottom of the domain was assumed to be fixed at saturation of 75% for Layer 3 [72].
All simulations in the current study lasted for 5 years, which was long enough to observe performance efficiency of the HGSHP system and shallow ground thermal behaviour due to the extraction and injection of geothermal energy.The maximum time-step was taken as 8640 s.
As the initial status of the ground at the testing site was unknown, using the same approach in the literature [56], a pre-modelling exercise with a duration of one year was conducted to obtain the initial distributions of pore-water pressure, thermal conductivity, and temperature profile in the ground.In the pre-modelling, the initial ground temperature profile was approximated using the analytical expression proposed by Hillel [34]: where T(y, t) is the soil temperature at time t (days) and depth y beneath the ground surface (m), T a is the constant ground temperature, which was taken as 10 • C, A 0 is the annual amplitude of the surface soil temperature, which was taken as 3.33 • C, t 0 is the lag time from arbitrary start date to the occurrence of the minimum soil temperature in a year, which was taken as 14 days, and d is a damping depth and calculated as follows: where D h is thermal diffusivity, which was taken as 0.01416 m 2 /d, and ω is 2π/365 (day − 1 ).The initial saturation of each soil layer was assumed to be 75% [72].
The above boundary conditions were applied to the same computational geometry in Fig. 7 without running the ground-loop pipes.Fig. 9 presents the last results of the pre-modelling, which were used as the initial pore-water pressure distribution (PW 1) and ground temperature profile (GT 1) for the 5-year-long simulations.

Parametric study
Various influential factors, including ground moisture transfer, building thermal load mode, buried depth of ground loops, and initial ground temperature profile were investigated through the following simulations listed in Table 8.
The effect of moisture transport in the ground on the performance of the HGSHP system was studied first.As shown in Fig. 9(a), in addition to the case of PW 1, another case with a fixed zero pore-water pressure distribution (PW 2) was assumed, i.e., the ground was saturated all the time.Correspondingly, the profiles of the soil thermal conductivity for the PW 1 and PW 2 scenarios are illustrated in Fig. 9(b).From the figure, the thermal conductivity of each soil in the PW 2 scenario (always in saturated status) was higher than that in the PW 1 scenario (in unsaturated status).
In terms of building thermal load, two cyclic heating/cooling modes (BL 1 and BL 2) were designed and shown in Fig. 10.For both modes, two sub-cases were considered and denoted by A and B. In BL 1, only the heating process was carried out by two periods, one was from January to April and the other was from October to December.The BL 1-A and BL 1-B represent 5% and 8% of the monthly total heating load of the site, respectively.In BL 2, the cooling process was included.Due to the lack of information regarding the site's actual cooling demand distribution with   time, it was assumed that the total cooling load was equal to the total heating load in BL 2 and that the cooling load was spread equally from May to September.The system is assumed to continuously operate under the designed thermal loads.
Considering the varied ground temperature with depth [30], two buried depths of ground-loop pipes were studied.In addition to the 3.0 m buried depth in Fig. 7 (H = 3.0 m), a different buried depth of 3.5 m (H = 3.5 m) was considered, and the other pipe configurations remained the same.
The undisturbed ground temperature plays significant role in the performance of HGSHP systems [45].To investigate the effect of a small rise in the undisturbed ground temperature on the HGSHP system's performance, along with GT 1, another ground temperature profile (GT 2) whose temperature was 1 ℃ higher than GT 1 at the same depth was designed, as shown in Fig. 9(c).It was assumed that the initial distributions of the pore-water pressure and soil thermal conductivity for the GT 2 scenario were the same as those for the GT 1 scenario.

Influence of ground moisture transfer 6.1.1. Soil saturation and thermal conductivity
Two scenarios (PW 1 and PW 2) were designed to study the influence of ground moisture transfer on the performance of the HGSHP system.For the PW 1 scenario, the moisture distribution and the related porewater pressure in the domain were transient and varying owing to the coupled effects between heat and moisture transfers under the prescribed boundary conditions.Whilst for the PW 2 scenario, a zero porewater pressure was fixed in the domain, which means that the ground was always saturated during the simulation, thus the thermal-hydraulic coupled problem would degenerate into a thermal problem in PW 2.
Variations in the saturation of soil at the observation point P (0.5, − 3.0, 0) for 5 years are illustrated in Fig. 11 70% and 75%.A steady annual cycle state after approximately two years can be observed owing to the annually prescribed climatic conditions and thermal load.In comparison, the value of saturation remained 100% all the time for the WP 2 scenario.Due to the periodic dynamics, the results in Year 5 (Day 1460 to Day 1825) were studied in detail and presented in Fig. 11(b).From the figure, the general trends of saturation at point P for the two sub-cases of PW 1 were similar.There were two stages that the saturation would increase as a result of the heating process and climatic conditions, one was from mid-January to mid-February and the other was from mid-May   to the end of August.Among the two sub-cases of PW 1, slight differences in saturation can be found from mid-February to early March and from October to December.There was BL 1-B > BL 1-A, as the increased thermal loads would lead to lower ground temperatures and hinder the vaporization of water.
Based on the soil saturation in Fig. 11(a) and assuming unchanged porosity and constant thermal conductivities of solid, water, and air, the variations in the soil thermal conductivity at point P can be calculated according to Equations ( 8) and ( 9), which are illustrated in Fig. 12(a).From the figure, the soil thermal conductivity exhibited the same periodic pattern as the soil saturation.
Taking the results in Year 5 as an example, as shown in Fig. 12(b), the soil thermal conductivity at point P in the PW 1 scenario varied between 0.716 W/m/K and 0.774 W/m/K, much lower than that in the PW 2 scenario (1.16 W/m/K), which is due to a lower soil saturation in PW 1 than PW 2.
Fig. 13(a) presents the development of ground temperature at point P for 5 years, which also demonstrated a steady annual cyclic state after about two years.The detailed variations in Year 5 under PW 1 and PW 2 are compared in Fig. 13(b), and there were similar trends in the ground temperature development for WP 1 and WP 2. The lowest and highest ground temperatures in the WP 1 and WP 2 scenarios both took place at the end of January and September, respectively.The ground temperature would noticeably increase with the non-heating process started from 1st May (Day 1580) due to the heat recovery of the ground.At the end of September (Day 1733), the highest ground temperature in the WP 1 scenario was less than 1 ℃ lower than the initial undisturbed value (12.15 ℃), while it was slightly higher by 0.1 ~ 0.2 ℃ in the WP 2 scenario.
Under the same thermal load (BL 1-A or B), the ground temperature in PW 2 was always higher than that in PW 1.Take BL 1-A for instance, the lowest and highest ground temperatures in PW 2 were 5.78 ℃ and 12.37 ℃, respectively, while in PW 1 they were 3.57 ℃ and 11. 51 ℃, respectively.The increment in the ground temperature was caused by the different soil thermal conductivities in WP 1 and WP 2 as shown in Fig. 12, which were closely related to the soil saturation as shown in Fig. 11 due to the moisture transfer in the ground.Moreover, the ground temperature in PW 2 stayed above 0 ℃ even under the thermal load BL 1-B, while for WP 1, only the lowest thermal load BL 1-A saw a ground temperature not below 0 ℃.With the higher thermal conductivity in WP 2 than that in WP 1, a higher thermal load can be supported by the current design of the HGSHP system.

Inlet and outlet fluid temperatures
Fluid temperature is an important indicator reflecting the performance of the HGSHP system.The inlet and outlet fluid temperatures during the 5-year simulations in the PW 1 and PW 2 scenarios are plotted in Fig. 14(a) and Fig. 15(a), respectively.Due to the annually periodic feature after about two years, the inlet and outlet fluid temperatures in Year 5 were taken out to analyse, which are shown in Fig. 14(b) and Fig. 15(b), respectively.Because the fluid temperature was directly dependent on the adjacent ground temperature, the development of the inlet and outlet fluid temperatures in Fig. 14(b) and Fig. 15(b) resembled the ground temperature as shown in Fig. 13(b).The lowest and highest fluid temperatures occurred at the end of January and September, respectively.Compared with the initial value (12 ℃), the highest fluid temperature would be around 1 ℃ lower in the PW 1 scenario but 0.3 ℃ higher in the PW 2 scenario.
It can be seen from Fig. 14(b) and Fig. 15(b) that the fluid temperature in the PW 2 scenario was always higher than that in the PW 1 scenario at the same thermal load.Taking BL 1-B as an example, the inlet and outlet fluid temperatures in the PW 2 scenario ranged from 0.91 ℃ to 12.32 ℃ and from 2.02 ℃ to 12.32 ℃, respectively.In comparison, the inlet and outlet fluid temperatures in the PW 1 scenario were from − 2.41 ℃ to 11.28 ℃ and from − 1.35 ℃ to 11.32 ℃, respectively.
In the simulations, pure water was used as the circulating fluid, the fluid temperature should be no less than 0 ℃ to avoid freezing.Therefore, without considering the ground moisture transfer or overestimating the ground saturation, such as in the WP 2 scenario, the current design could work well under the thermal loads of BL 1-A and B; in contrast, for the realistic condition that considered the ground moisture transfer, such as in the WP 1 scenario, the current design of the HGSHP system would only allow the thermal load of BL 1-A.

Heat pump COP
The heat pump COP is a direct indicator showing the performance efficiency of the HGSHP system.Variations in the heat pump COP for 5 years and in Year 5 in the PW 1 and PW 2 scenarios are illustrated in Fig. 16(a) and Fig. 16(b), respectively.Taking the results in Year 5 to analyse, as shown in Fig. 16(b), for both PW 1 and PW 2, the heat pump COP was higher than 3.6 and lower than 5.0.Additionally, the greatest COP all took place on Day 1733 (the end of September), indicating the non-heating period from May to September was beneficial to the performance of the HGSHP system because the ground temperature would recover with time.However, the heat pump COP was higher in WP 2 than in WP 1 under two thermal loads due to the higher soil saturation and thermal conductivity in WP 2 than in WP 1.
Regarding the effects of building thermal load mode, buried depth of ground loops, and initial ground temperature profile, the results in Year 5 will be illustrated and analysed since an annual periodic state after about two years was found in each scenario.

Influence of building thermal load mode
The ground temperature at point P (0.52, − 3.0, 0), inlet and outlet fluid temperatures, and heat pump COP under two building load modes (BL 1 and BL 2) in Year 5 (Day 1460 to Day 1825) are shown in Fig. 17.In the simulations, the buried depth of pipes was 3.0 m (H = 3.0 m), and the initial pore-water pressure distribution PW 1 and the initial ground temperature profile GT 1 were adopted.
As mentioned earlier, unlike BL 1 that was only a heating process, BL 2 included both heating and cooling processes.From Fig. 17   would exceed the undisturbed ground temperature on Day 1587 and continue to increase till the cooling process ended.Taking the case of BL 2-B as an example, its highest ground temperature at point P reached 37.77 ℃, much higher than the location's undisturbed ground temperature of 12.14 ℃ on Day 1733.However, as the heating process initiated after Day 1733, a more rapid decline in the ground temperature can be found in the BL 2 scenario compared with the BL 1 scenario due to the higher thermal gradient.Among the two sub-cases of BL 2, although the amount of extracted and injected heat was smaller for the case of BL 2-A, the ground temperature of BL 2-A was higher than BL 2-B from Day 1460 to Day 1580 (the first heating period from January to April) and from Day 1764 to Day 1825 (October to December).This implied that the injection of large amounts of heat into the shallow ground via horizontal ground heat exchangers cannot offset the reduction in ground temperature due to heat extraction, since most of the injected heat would be dissipated.
Fig. 17(b) and Fig. 17(c) illustrate the inlet and outlet fluid temperatures, respectively.The fluid temperature showed a similar development trend as the ground temperature in Fig. 17(a).Owing to the cooling process, the fluid temperature in BL 2 was higher than that in BL 1 in each sub-case.For instance, the lowest inlet and outlet temperatures for BL 2-A occurred at the end of January (Day 1491) with the values of 3.35 ℃ and 3.97 ℃, respectively.In comparison, for BL 1-A, the lowest inlet and outlet temperatures took place on the same day but were 2.38 ℃ and 3.08 ℃, respectively.As the cooling process started in May (Day 1580), the fluid temperature soon exceeded its initial value in a couple of days and peaked at the end of September (Day 1733).For example, the highest inlet fluid temperatures for BL 2-A and BL 2-B were 29.87 ℃ and 41.80 ℃, respectively.However, when the heating process carried out, the fluid temperature decreased rapidly.More explicitly, as the heating process carried on, the fluid temperatures in BL 2-B would drop below the freezing point of water from January to April, which demonstrated the failure of the HGSHP system.
The COP values under different building thermal loads are shown in Fig. 17(d), which were greater than 3.25.From January to April (the first heating period of the year), small differences in the COP values can be found between BL 1 and BL 2. From October to December (the second heating period of the year), the COP values in BL 2 were greater than that in BL 1 due to the heat injection from May to September (Day 1580 to Day 1733).In addition, BL 2-A demonstrated higher COP values than BL 2-B all year round except October (Day 1733 to Day 1764) although the lowest amount of heat was dumped.
The above results indicated that the cooling process could improve the performance of the HGSHP system to a certain extent during the heating periods, since most of the injected heat would be dissipated.

Influence of buried depth of ground loops
Fig. 18 shows the ground temperatures at point P (0.52, − 3.0, 0) and (0.52, − 3.5, 0), inlet and outlet fluid temperatures, and heat pump COP in Year 5 when the pipes were buried at two different depths (H = 3.0 m and H = 3.5 m).In these two scenarios, the initial pore-water pressure distribution PW 1 and the initial ground temperature profile GT 1 were employed, and the HGSHP system operated under the building thermal load mode BL 1.
The ground temperatures at the locations of (0.52, − 3.0, 0) and (0.52, − 3.5, 0) are compared in Fig. 18(a), in which the varied undisturbed ground temperatures at the depths of 3.0 m and 3.5 m are plotted as well.Based on the initial ground temperature profile (GT 1) in Fig. 9 (c), the initial ground temperature at point P in the H = 3.5 m scenario was 11.93 ℃, which was 0.12 ℃ higher than that in the H = 3.0 m scenario (11.81 ℃).The average undisturbed ground temperature in the H = 3.5 m scenario was 11.87 ℃ (11.59 ℃ to 12.14 ℃), whereas it was 11.73 ℃ in the H = 3.5 m scenario (11.17 ℃ to 12.28 ℃).When geothermal energy was extracted from the ground, owing to the higher initial and average ground temperatures at H = 3.5 m than H = 3.0 m, the ground temperature at (0.52, − 3.5, 0) was greater than that at (0.52, − 3.0, 0) under each sub-case.
The inlet and outlet fluid temperatures are shown in Fig. 18(b) and Fig. 18(c), respectively.Due to the ground temperature differences, the fluid temperature in the H = 3.5 m scenario was higher than that in the H = 3.0 m scenario.Taking the thermal load of BL 1-B as an example, the lowest inlet and outlet fluid temperatures in the H = 3.0 m scenario both fell below 0 ℃, while they were higher 0 ℃ in the H = 3.5 m scenario.The results indicated the HGSHP system can operate properly at a higher thermal load when the ground loops were buried at 3.5 m instead of 3.0 m.Fig. 18(d) shows the heat pump COP when ground loops were buried at two depths.Under the same thermal load mode, a higher value of heat pump COP was obtained at H = 3.5 m compared with H = 3.0 m.For instance, under the thermal load of BL 1-A, the heat pump COP varied from 4.86 to 4.24 at H = 3.5 m and from 4.83 to 4.00 at H = 3.0 m.
Since the shallow ground temperature is significantly affected by climatic conditions, the buried depth plays an important role in the proper operation of the HGSHP system.Accurate measurement or assessment of the shallow ground temperature should be conducted before the design of the HGSHP system.

Influence of initial ground temperature profile
Under two initial ground temperature profiles (GT 1 and GT 2), the ground temperature at point P (0.52, − 3.0, 0), the inlet and outlet fluid temperature, and the COP of the heat pump in Year 5 are presented in Fig. 19.In these two scenarios, the pipes were buried at depth of 3.0 m with the pore-water pressure distribution PW 1, and the HGSHP system worked under the building thermal load BL 1.
Fig. 19(a) shows variations in the ground temperature at point P using GT 1 and GT 2, and the same development trends can be identified.As the initial ground temperature in GT 2 was 1 ℃ higher than that in GT 1 at the same depth, the ground temperature of GT 2 was at least 0.75 ℃, and 0.71 ℃ higher than those of GT 1 under the thermal load BL 1-A and BL 1-B, respectively.Moreover, for the scenario of GT 1, the ground temperature at point P can only stay above 0 ℃ under BL 1-A, while for GT 2, a higher building load BL 1-B saw a ground temperature higher than 0 ℃ (the lowest value was 0.14 ℃ took place at the end of January).
Fig. 19(b) and Fig. 19(c) present the inlet and outlet fluid temperatures, respectively.In agreement with the ground temperature, the inlet and outlet fluid temperatures of GT 2 were higher than those of GT 1.
Developments of COP for GT 1 and GT 2 are shown in Fig. 19(d).Due to the higher outlet fluid temperature of GT 2, the value of COP for GT 2 was at least 0.05 higher than that for GT 1, implying better operational performance of the HGSHP system under GT 2 than GT 1.
The results showed that a 1 ℃ rise in the initial ground temperature can improve the performance efficiency of the HGSHP system.The HGSHP system may operate sustainably under a higher building load, such as BL 1-B if a different fluid with a lower freezing point was utilized

Future work
When the HGSHP system is subjected to high heating load, the ground temperature around the pipes may drop below the freezing point of water (0 ℃), and the liquid water in soil voids will form ice.In order to ensure the service life of the buried pipes and the efficiency and safety of the system, ground freezing should be limited or avoided because the pipes can be squeezed by the volume expansion in freezing soil [60].Ice in soil voids will thaw as ground temperature recovers between heat extraction periods or due to heat injection.The phase change from liquid water to ice and vice versa changes the physical properties of water substantially [1].In future work, the effects of freeze-thaw of the ground will be considered in the coupled TH model.

Conclusions
This paper presents a 3D coupled TH model with realistic boundary conditions to study the performance of a HGSHP system and the thermal behaviour of the shallow ground in response to the extraction and injection of geothermal energy.The ground temperature profile was obtained by the coupled TH model based on the real climatic conditions with the consideration of soil moisture variation.The ground loops and the heat pump were coupled by the heat pump COP.The proposed model was validated against experimental results and an analytical model obtained from literatures.
Long-term numerical simulations were performed to investigate the performance efficiency of a HGSHP system to meet the thermal demand of a potential testing site for 5th generation district heating and cooling network development.Performance of the HGSHP system and the ground thermal behaviour were predicted and evaluated by accounting the whole system including building heating and cooling load, heat pump, ground heat exchangers, and the ground.
Influencing factors, including ground moisture transfer, building thermal load mode, buried depth of ground loops, and initial ground temperature profile were studied to assess performance efficiency of the HGSHP system and evolution of ground thermal behaviour in response to geothermal energy extraction or rejection into ground.Based on the simulations, the following conclusions can be drawn: • The performance of the HGSHP system is significantly influenced by moisture transfer in the ground.Overlooking the ground moisture transfer or hyperbolizing the ground saturation would overestimate the load-carrying capacity of the HGSHP system.• The cooling process can improve the performance of the HGSHP system relative to the heating-only mode, whilst it is not appropriate for long-term thermal storage of the HGSHP system due to fast heat dissipation.• Buried depth of ground heat exchangers has a considerable impact on the performance of the HGSHP system.Heat imbalance tends to be serious if the ground heat exchangers are not buried deep enough.• Given the heat carrier fluid is pure water, neither the ground temperature nor the inlet fluid temperature should fall below 0 ℃ to prevent freezing of both the ground and the circulating fluid.For the testing site, 5% of the monthly total heat demand could be met by the designed HGSHP system, consisting of 200 U-shaped ground loops buried at the depth of 3 m.The maximum heat demand for the site in

06 Fig. 5 .
Fig. 5. Comparison of the outlet fluid temperature predicted by the presented model and reported in Hu [35].

0075 Fig. 6 .
Fig. 6.Comparison of the outlet fluid temperature predicted by the presented model and reported in Li et al. [44].

Fig. 7 .
Fig. 7. Representative model simulation domain with initial and boundary conditions as well as 3D discretized finite element mesh.

Fig. 8 .
Fig. 8. Atmospheric data obtained from Met Office [46]: (a) ambient air temperature T a , (b) shortwave solar radiation H SW , (c) air relative humidity h a , (d) rainfall P, and (e) wind speed..u y (a).As shown in the figure, the saturation value at point P in PW 1 generally fluctuated between

Fig. 10 .
Fig. 10.Building thermal load per U-shaped ground-loop: (a) BL 1 -heating process and (b) BL 2 -heating and cooling processes (the value of heating load is positive, and the value of cooling load is negative).

Fig. 14 .
Fig. 14.Variations in inlet fluid temperature: (a) for 5 years and (b) in Year 5 under two pore-water pressure distributions PW 1 and PW 2.

Fig. 15 .
Fig. 15.Variations in outlet fluid temperature: (a) for 5 years and (b) in Year 5 under two pore-water pressure distributions PW 1 and PW 2.

Fig. 16 .
Fig. 16.Variations in heat pump COP: (a) for 5 years and (b) in Year 5 under two pore-water pressure distributions PW 1 and PW 2.

Table 5
Physical parameters of soil layers.
W.Gao et al.

Table 6
Parameters of ground-loop pipes.

Table 7
Parameters of fluid in pipes.

Table 8
Summary of simulations conducted in this study.