Improving computational efficiency of numerical modelling of horizontal ground source heat pump systems for accommodating complex and realistic atmospheric processes

Modelling horizontal ground loops for a horizontal ground source heat pump (HGSHP) system is complex and computationally expensive. The computation precision is highly reliant on the prescription of an undisturbed ground temperature in the unsaturated ground as well as realistic and accurate atmospheric processes at the ground surface boundary. Conventionally, modelling of such a system would include direct application of the atmospheric processes at the soil-atmosphere boundary and solve it in a single-stage approach. However, low efficiency is found for large spatial domain and long-term transient problems as the boundary processes need to be solved and expressed in terms of primary model variables at each simulation time-step. This paper proposes an equivalent two-stage modelling approach, for the first time, based on an advanced coupled thermal-hydraulic (TH) model to improve computation efficiency while maintaining adequate accuracy. In this approach, firstly, the model is solved for an intact ground that is imposed by complex atmospheric processes, e.g., rainfall, solar radiation, humidity, evaporation, etc. at the soil-atmosphere boundary, and the spatial and temporal variations of the primary model variables are recorded. Afterwards, the recorded data are incorporated in the simulator, as model inputs, for the same ground including a HGSHP system. Predicted results from both 2D and 3D simulations show that the ground temperatures calculated by the proposed two-stage approach are in good agreement with that of the traditional single-stage approach. However, the two-stage approach is computationally robust. For the presented 2D and 3D simulations, it required only 32% and 37% of the time of the single-stage approach, respectively, while maintaining great accuracy. This demonstrates the utility of the proposed two-stage approach for modelling complex scenarios of realistic HGSHP systems installed in a large spatial domain and for long-term operation.


Introduction
In recent years, the ground source heat pump (GSHP) system has received increasing attention, which circulates a heat carrier fluid, usually in closed-loops, to extract geothermal energy to heat buildings or inject extra heat underground to cool buildings.Horizontal GSHP or HGSHP systems that use horizontal ground-loops as ground heat exchangers (GHEs) are very common due to the low initial installation cost and ease of balance between thermal efficiency and expense (Li et al., 2017;Kim et al., 2018).Numerical modelling is widely used in long-term investigations of HGSHP system performance and ground thermal behaviour, as well as in the design of HGSHP systems, due to its flexibility for analysing various scenarios, and direct demonstration of predicted results or performance.However, numerical modelling of HGSHP systems accounting realistic ground processes and, more importantly, implementing true atmospheric processes at the soil-atmosphere boundary are computationally expensive and challenging.Improvement of modelling strategies and/or implementation methodologies is therefore essential to address the issues.
Several studies that focused on modelling HGSHP systems are available in literatures.For example, Esen et al. (2007) built a numerical model to simulate the ground temperature distribution in the vicinity of pipes with the operation of a HGSHP system and compared the simulated results against the monitored data.They concluded that the heating load of the HGSHP depends on the ground temperature distribution around the pipes.Dasare and Saha (2015) developed a numerical model to evaluate the thermal performance of various types of horizontal GHEs, viz.linear, helical, and slinky, and the influences of buried depth, fluid flow rate, and soil thermal conductivity were investigated.They found that the thermal conductivity of soil is the most important parameter in the heat transfer process and the helical geometry is the best performing amongst the other geometries.Li et al. (2017) established a numerical model to study the operation characteristics of the horizontal spiral-coil GSHP system, including the inclusion of heat pump, subsurface factors, daily variations of load, and operation models.They found the importance of the heat pump Coefficient of Performance (COP) on the fluid temperature, soil thermal conductivity and pipe spacing are main factors to influence the system's performance, and continuous operation shows the best performance.Based on the ANSYS Fluent software, Pu et al. (2018) evaluated the performance of various pipe arrangements (in-line and staggered) and investigated the effects of the relative displacement of staggered pipes.Moreover, the effects of bending number, pipe spacing, and buried depth were studied.The results showed that the staggered pipes outperformed the in-line arrangement when the relative offset displacement was greater than 1/3, and the critical pipe spacing for doubled-layered horizontal GHEs was obtained.Sedaghat et al. (2020) developed a Computational Fluid Dynamic (CFD) model to evaluate the long-term performance of a novel HGSHP system for space cooling in a hot climate and investigated the influences of pipe configurations on the system performance.They reported that not only the average annual COP but also the maximum cooling load can be improved by the novel system.Using the FEFLOW software, Bina et al. (2020) conducted a sensitivity analysis to observe the influences of various design parameters on the performance of a horizontal GHE drilled by the horizontal directional drilling technology.It showed that the heat exchange rates were greatly enhanced in large diameter system but declined in longer length system, and faster groundwater velocities in perpendicular direction can improve the system performance.Tang and Nowamooz (2020) proposed a numerical framework considering the atmosphere-soil-HGHE interaction, and COMSOL software was used to evaluate the outlet temperatures of a slinky-type HGHE installed in a multi-layered soil field by considering the local metrological and geological conditions.The results showed that the increase of the installation depth increased the outlet temperatures, and non-consideration of the atmosphere-soil interaction overestimated the annual fluid outlet temperature in the heating scenario.Also employing COMSOL software, Shi et al. (2022) numerically analysed influential factors, including injection mass flow rate, operation modes, heat storage, subsurface water flow, soil thermal conductivity, installation depth and atmospheric conditions (precipitation and air temperature), on the HGHE performance.The results showed subsurface water flow, injection mass flow rate and climatic conditions should be priorities for the optimal design of HGHE.Xu et al. (2022) built a 3D CFD model with single layer ground heat exchanger to simulate the sandbox experiment to investigate the extent of soil thermal imbalance and its effects on performance of HGSHP system.It was found that optimized intermittent mode and increasing soil moisture content are effective solutions to eliminating thermal accumulation around buried pipes and improving the performance of HGSHP systems.
In addition, numerical modelling is adopted by researchers to validate the design methods for HGSHP systems.For instance, Kim et al. (2018) proposed a design equation for the design length of a horizontal spiral-coil GHE, and the CFD program COMSOL was utilized to verify the proposed design method.
From various perspectives, researchers have proposed different methods and techniques to improve the computation efficiency of modelling GSHP systems, especially for large systems coupled with vertical boreholes.For instance, Kim et al. (2010) presented a numerical model for short-time transient response of heat transfer around boreholes, using a state model size reduction technique for decreasing the calculation time.Cullin and Spitler (2011) developed a methodology for determining the duration of the shorter timestep and the magnitude of the corresponding load, which increases the performance of heat exchanger design tools that incorporate a "hybrid timestep" approach.Brunetti et al. (2017) proposed a pseudo-three-dimensional (3D) model, which combines a one-dimensional (1D) description of the heat transfer in the buried tubes of the exchanger with a two-dimensional (2D) description of the heat transfer and water flow in the surrounding subsurface soil, to reduce the computational cost for the numerical analysis and interpretation of thermal response tests (TRTs).Fang et al. (2018) developed a software package based on the finite difference method for thermal analysis of deep borehole heat exchangers (DBHEs), and the computation efficiency is increased by using an algorithm for the direct solution of resulted algebraic equation set.Zhang et al. (2021) established a heat transfer model for a deep borehole cluster and proposed a new dimension reduction algorithm based on the linear superposition principle to improve computation efficiency in solving complex 3D heat transfer problems.Moreover, Huang et al. (2022) developed two approximation-assisted reduced-order modelling methods for a phase change thermal storage device, including a pure black-box model and a grey-box model based on the Number of Transfer Units (or effectiveness-NTU) approach, and much lower computational time than finite-volume models was required.
As the ground is worked as a thermal source/sink, the performance and thermal efficiency of HGSHP systems are significantly influenced by the undisturbed ground temperature (Radioti et al., 2017;Jeong et al., 2017).In the ground, the undisturbed ground temperature within a certain depth range from the ground surface fluctuates with time because of the influences of the local climatic conditions and the ground structure conditions.This depth is from 0 m to about 8 m to 20 m (Bryś et al., 2018).The ground-loops of the HGHSP system are usually buried at 1-3 m underground (Li et al., 2017;Kim et al., 2018).Therefore, the undisturbed ground temperature that changes temporally and spatially is especially of importance for the performance of HGSHP systems.
In the numerical modelling of a HSGHP system, prescribing the undisturbed ground temperature is one of the key factors that should be realistic.In numerical modelling literature (Sanaye and Niroomand, 2010;Dasare and Saha, 2015;Li et al., 2017;Kim et al., 2018;Pu et al., 2018), a sinusoidal equation to explicitly characterize the ground temperature as a function of depth and time was used.It should be noted that this analytical expression was derived by Kusuda and Achenbach (1965) for a semi-infinite solid considering heat conduction below the ground surface owing to the ambient temperature.Although this method is beneficial to increasing the numerical modelling computation efficiency, the undisturbed ground temperature is roughly dictated during numerical modelling, increasing the cost of computation accuracy.In contrast, in the numerical model for a HGSHP system proposed by Kayaci and Demir (2018), the transient ground temperature profile was obtained by the real climatic conditions on the ground surface boundary.However, this model was built only on the heat transfer and ignored the moisture transfer in soils that are often in the unsaturated status.Moreover, soil thermal properties are not constant but significantly influenced by its water content (saturation), which has been revealed by experiments.For example, the thermal conductivity of sandy soil can be reduced from 2.65 W/m/K in saturated condition to 0.9 W/m/K in dry soil (Akrouch et al., 2015).The thermal conductivity of Coode Island Silt, which is a silty clay, was decreased from 1.34 W/m/K to 0.34 W/m/K when its saturation decreased from 100% to 0% ( Barry-Macaulay et al., 2013).In addition, to a better knowledge of the relationship between thermal conductivity and water content, required for understanding the behaviour of HGSHP systems and enhancing the system performance, based on the Improving Thermal Efficiency of horizontal ground heat exchangers (ITER) Project, Di Sipio and Bertermann (2017Bertermann ( , E. 2018aBertermann ( , E. 2018b) conducted field tests of horizontal helix earth collectors in the same climatic conditions and under the same thermal stress for five different soil mixture.It revealed that the change of soil moisture content is the main parameter influencing the thermal properties of the soils, and compared with coarse sand, loamy sand and W. Gao et al. fine sand can provide better performance for HGSHP systems due to the high capacity of hosting water, which significantly improves the thermal conductivity and the latent heat storage capacity.Gonzalez et al. (2012) also conducted field observations for a HGSHP system in the south of the UK and summarised the effect of heat extraction on the soil physical environment.Soil temperatures and soil moisture content measurements showed that the system performance was affected by the heat and water transport in the soil.Overall, such a numerical model that can obtain a realistic ground temperature profile that changes transiently, as well as consider the heat and moisture transfer in unsaturated soils and dynamic soil thermal properties, is still lacking in the literature.
To improve numerical modelling performances, a new modelling procedure is proposed in this paper.An advanced coupled model that studies both thermal and hydraulic (TH) responses of a ground that facilitate horizontal ground heat loops and subjected to realistic atmospheric and ground conditions is developed.The primary model variables are temperature (T) and porewater pressure (u).The atmospheric conditions include solar radiation, rainfall, humidity, air temperature, and wind velocity, and the ground conditions include ground water and temperature flow.In a traditional approach, modelling of HGSHP systems would include direct application of the atmospheric processes at the soil-atmosphere boundary and solve it using a numerical simulator in a single-stage approach.However, for large spatial domain and longterm transient problems, model simulations are found to be computationally expensive, as the boundary processes need to be solved and estimated in terms of the primary model variables, e.g., T and u at each simulation time-step.This issue has been addressed in this study following an equivalent two-stage simulation approach.Firstly, the model is solved for an intact ground imposing the atmospheric boundary conditions, and the spatial and temporal variations of the primary variables on the ground surface are recorded.Afterwards, the model variables are incorporated in the simulator for the same ground using a HGSHP system.Accuracy of this equivalent two-stage approach and improvement of computational times are investigated comprehensively by conducting both 2D and 3D model simulations.

Coupled thermal-hydraulic model for unsaturated soil
Moisture flow within unsaturated soil is described as a two-phase process, consisting of water and vapour flows.The general equation for the moisture flow can be described as: in which θ 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 water, ρ v is the density of vapour, v l is the velocity of water, and v v is the velocity of vapour.Conduction and convection are the two heat transfer mechanisms considered within unsaturated soil.The temporal derivative of the heat content is equal to the spatial derivative of the heat flux, which can be expressed as follows: in which L is the latent heat of vaporisation, T is the temperature, ϕ is the soil porosity, T r is the reference temperature, 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, ρ s is the density of the solid, and λ T is the thermal conductivity.
The water velocity is described by Darcy's law, and the vapour velocity is obtained via the equation proposed by Philip and de Vries (1957): where u is the pore-water pressure, γ l is the unit weight of water, y is the elevation, K l is the unsaturated hydraulic conductivity, 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.
The Brooks and Corey (1964) Model is adopted to model the unsaturated hydraulic conductivity, and there is: where K ls is the saturated hydraulic conductivity, θ ls is the saturated water content, and η is the shape parameter.The Van Genuchten (1980) Model is used to characterize the Soil Water Characteristic Curve (SWCC) of soils, which is expressed by: Fig. 1.Schematic of horizontal ground source heat pump system and climatic conditions on the ground surface boundary.
W. Gao et al. ( where h is the pressure head, h g is the scale parameter, and n is the shape parameter. Based on the unsaturated soil's composition (Thomas and Rees, 2009), its thermal conductivity can be obtained: 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:

Ground surface boundary
Fig. 1 shows the HGSHP system installed in the shallow ground.Under the ground, multiple U-shaped ground-loops that buried at the same depth form the ground heat exchangers.On the ground surface, namely, at the soil and atmosphere interface, energy and moisture exchanges occur constantly under the influence of climatic conditions, including solar radiation, rainfall and clouds, air temperature and humidity, and wind.The climatic conditions are assumed to be distributed uniformly on the surface of the ground.As a result, the temperature of the shallow ground exhibits seasonal variations with depth, which will influence the heat exchange between the horizontal ground loops and the adjacent ground, hence affecting the performance of the HGSHP system.
The energy balance equation can be developed at the ground surface boundary (Mihalakakou et al., 1997;Bryś et al., 2018), which is given as follows: where H E 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.
For the absorbed shortwave radiation flux at the ground surface H Absorbed SW , it can be expressed by the following equation (Deardorff, 1978): where ε SW is the shortwave reflection factor associated with the ground surface type, and H SW is the shortwave solar radiation.
The longwave radiation flux absorbed at the ground surface H Absorbed LW is given as (Imberger and Patterson, 1981;Lewis et al., 2004): in which σ is the Stefan-Boltzmann constant (5.67 × 10 − 8 W/m 2 /K 4 ), C cloud is the fractional cloud cover coefficient (0 for clear sky and 1 for totally overcast), T air 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 can be obtained by the following equation: Based on the Stefan-Boltzmann's law (Woodward et al., 2001), the longwave radiation flux emitted from the ground surface H Emitted LW is calculated: where ε LW is the long-wave emissivity of the ground.
For the sensible heat radiation flux at the ground surface H SEN , the following relation is used (Deardorff, 1978): 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.
The latent heat radiation flux at the ground surface due to evaporation is obtained by the equation below (Deardorff, 1978): where L is the latent heat of vaporization, and E is the evaporation flux (saturated state), which can be calculated as follows (Sverdrup, 1946;Deardorff, 1978): in which q is the specific humidity of the soil at the ground surface, and q air is the specific humidity of air.As the ground surface would enter the unsaturated state, the following modification is made to Eq. ( 16) and there is (Barton, 1979;Wilson et al., 1997): where E A is the actual evaporation flux (unsaturated state), h ground is the relative humidity of the ground surface, and h air the relative humidity of air at the ground surface.
Assuming the moisture flux at the ground surface is a hydrological process, the balance equation in terms of moisture can be presented as follows (Fredlund et al., 2011): where Q M is the net moisture flux at the ground surface, P is the rainfall, E A is the evaporation flux, and R O is the run-off.By collectively observing Eqs. ( 9)-( 18), five climatic variables are needed to model the specific ground surface boundary in terms of energy and moisture transfer of the coupled thermal-hydraulic (TH) model, and they are ambient air temperature T air , shortwave solar radiation H SW , air relative humidity h air , wind speed u y , and rainfall P.These climatic variables are commonly monitored in the field and representative values can be obtained for most regions globally from metrological data.

Ground-loop boundary
When the HGSHP system is operating, the circulating fluid exchanges heat with the adjacent ground through the pipe wall, thus heat flux is applied on the ground-loop boundary.Considering the fluid is normally turbulent in the pipes, and the pipes are in good contact with the ground, the pipe wall thickness is negligible.In addition, it is assumed that the distance between pipes is big enough to avoid thermal interference between them.
In 2D modelling, the heat flux applied at the ground-loop boundary can be assumed to be uniform for simplification, i.e., the heat flux is equal on the perimeter of pipes.The heat flux per unit area of a groundloop Q A can be expressed as: where Q GL is the thermal load of a ground-loop, and L GL is the total length of a ground-loop.In 3D modelling, considering the thermal gradient of the fluid along its circulating direction, the heat flux on the ground-loop is nonuniform.The fluid temperature profile should be predicted first, and then the heat flux at the ground-loop boundary can be calculated based on the locally transient fluid and ground temperatures.To obtain the fluid temperature profile, the ground-loop is discretized into a series of control volumes with a length of dL, and the fluid temperature in each control volume T f is assumed to be constant.The heat flux per area of a ground-loop Q A can be given by the following equation: where λ f is the thermal conductivity of the fluid, T g is the ground temperature adjacent to the control volume, and R is the pipe radius.

Numerical solution and equivalent two-stage approach to represent climatic variables
The presented coupled TH model has been developed within the computer code COMPASS (Code of Modelling Partially Saturated Soils), which is developed at the Geoenvironmental Research Centre, Cardiff University (Thomas and He, 1995).In the model, the governing equations are expressed in terms of the primary variables (pore-water pressure u and temperature T) as follows: where C and K terms represent storage and flux, respectively.For the detailed expression of each term, please refer to the literature (Gao et al., 2022).Within the framework of COMPASS, the above equations are spatially discretised using Galerkin finite-element method (Zienkiewicz et al., 2005), 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 (Douglas and Jones, 1963) to obtain the converged primary variables.Both the ground surface boundary and the ground-loop boundary are implemented in COMPASS.The coupled TH model with the sophisticated boundary conditions has been verified and validated against experimental results and analytical model (Hepburn, 2013;Gao et al., 2022).
During the computation of the numerical model, the inputs of climatic conditions are in fact converted into pore-water pressures and temperatures at the nodes of the model domain surface.However, this conversion process can be time-consuming as the domain size of a HGSHP system is often large and the operation time is usually long.Rather than executing large numerical simulations with real atmospheric processes at each time-step, i.e., the traditional single-stage approach, it is proposed that equivalent pore-water pressures (u) and temperature (T) at surface nodes be obtained by modelling an intact ground with the real climatic conditions first.The spatial and temporal variations of the primary variables on the ground surface are recorded.After that, the equivalent u and T are substituted for the atmospheric conditions in a large modelling of a HGSHP system.This approach is hereinafter referred to as the equivalent two-stage approach.

Model domain
In this study, both 2D and 3D simulations were carried out using the traditional single-stage approach and equivalent two-stage approach to compare the ground temperatures with the operation of a HGSHP system.

Model parameters and thermal load
The physical parameters of three soil layers are listed in Table 1, and the parameter values were determined from data in literatures (Leij et al., 1996;Busby, 2015;Song and Hong, 2020;Parkes et al., 2021).The physical parameters of fluid in pipes are provided for 3D simulations and listed in Table 2. Pure water was chosen as the heat carrier fluid due to its non-toxicity, low expense, and good thermal properties (Bartolini et al., 2020;Soltani et al., 2021).The initial fluid temperature was

Table 1
Physical parameters of three soil layers (Gao et al., 2022).determined based on the cold running water in the UK and the undisturbed ground temperature.In the simulations, the fluid was assumed to enter from the right pipe and exit from the left pipe.Moreover, the flow rate was assumed to the same in each pipe with a value of 2.5 m 3 /h.Thus, a turbulent flow with a Reynolds number of approximately 12,400 was created, followed by a sufficient convective heat transfer between the fluid and the adjacent ground.
The thermal load of the ground-loop used in 3D simulations is shown in Fig. 3a.As can be seen in the figure, the heating process included two periods in a year, one was from January to April, and the other was from October to December.The corresponding heat flux used in 2D simulations, as shown in Fig. 3b, were computed based on Eq. ( 19) and the pipe configuration.

Initial conditions and boundary conditions
As mentioned earlier, five climatic variables are needed to model the ground surface boundary.The climatic variables in the 2D and 3D simulations are presented in Fig. 4, which were monitored at a meteorological station at Church Lawford, UK in 2019 (Met-Office, 2021).The other ground surface boundary-related parameters were obtained from the literature (Van Wijk, 1966;Calvet et al., 1999;Mihalakakou et al., 1997;Hepburn, 2013;Staniec and Nowak, 2016), as listed in Table 3.
For the 2D and 3D domains as shown in Fig. 2, the left and right sides were symmetrical.At the domain bottom, a fixed temperature boundary (285.15K) was prescribed to consider the undisturbed ground temperature and the pore-water pressure was assumed to be fixed at saturation of 0.75.Fig. 5 shows the initial temperature and pore-wate pressure along the buried depth.
It should be pointed out different maximum time-steps were adopted  for the two approaches.The maximum time-step for simulations using the traditional single-stage approach was 1728 s for convergency, while for the simulations using the equivalent two-stage approach, a higher maximum time-step (5000 s) can be prescribed and the convergency can still be ensured.

Temperature and pore-water pressure equivalent to the climatic conditions
By performing a simulation with a 2.0 m × 2.0 m × 4.0 m domain (W × L × D) without operating the HGSHP system, the temperature and pore-water pressure at the nodes of the model domain surface, which are equivalent to the climatic variables in Fig. 4, are illustrated in Fig. 6.From the figure, the temperature of ground surface varied from 275. 3 K to 291.1 K in a year, and the pore-water pressure of ground surface ranged from -1.77 × 10 4 Pa to -2.24 × 10 7 Pa in a year.Ground surface temperature and pore-water pressure both changed month by month as a result of the monthly average climatic conditions as shown in Fig. 4. The period of the highest ground surface temperature coincided with the period of the lowest ground surface pore-water pressure.Furthermore, the stage of rising temperature corresponded to the stage of decreasing pore-water pressure, and vice versa.

2D modelling
Using the monitored climatic conditions (seen in Fig. 4) and the equivalent pore-water pressure and temperature at surface nodes (seen in Fig. 6), respectively, 5-year-long 2D simulations were conducted.The simulated ground temperatures at P1, P2, P3, and P4 are compared in Fig. 7. From the figure, it can be seen that the ground temperature exhibited a periodic cycle from the second year regardless of the approach to set the ground surface boundary.Based on the results of Year 3 at P1, P2, and P3, the deeper the location, the smaller the ground temperature amplitude change.At P4, which was next to the right pipe where cold fluid entered, the ground temperature showed greater temperature variations than that of P3, which had the same depth as P4.
Fig. 7 shows that the ground temperatures calculated by the twostage approach were in good agreement with the results obtained by the single-stage approach.The differences in the ground temperatures between the two approaches decreased with the increase of the buried depth.For example, at P1 (Fig. 7a), the ground temperature varied from 277.3 K in winter to 285.9 K in autumn by the single-stage approach, while it was from 277.9 K in winter to 287.1 K in autumn by the twostage approach, creating a maximum difference of 1.2 • C (K).At P3 (Fig. 7c), the maximum difference in the ground temperature between the single-stage approach (277.6K to 284.5 K) and the two-stage approach (277.6K to 284.8 K) reduced to 0.3 • C (K). Especially, at the W. Gao et al. location of P4 (Fig. 7d), the ground temperatures obtained by the two approaches basically coincided.

3D modelling
3D simulations that last for 5 years were performed by the singlestage approach and the two-stage approach, respectively.Fig. 8 illustrates the ground temperatures at M1, M2, M3, and M4 in the 3D domain.From the figure, the ground temperatures at various depths obtained by both approaches reached a steady annual cyclic state from the second year due to annually repeated boundary conditions.The deeper the location, the less the ground temperature amplitude altered, according to the results of Year 3 at M1, M2, and M3.In addition, owing to rapid geothermal energy extraction from the ground adjacent to the pipe, more distinct variations in the ground temperature can be observed at M4.
Upon the collective inspection of Fig. 8, it can be seen that the twostage approach's ground temperatures closely matched the results simulated by the single-stage approach, with a maximum difference of 0.2-0.6 • C (K) at various depths.For instance, the ground temperature at M1 (Fig. 8a) by the single-stage approach ranged from 278.3 K in winter to 287.3 K in autumn, while from 277.8 K in winter to 287.2 K in autumn by the two-stage approach.The ground temperature at the location of M3 (Fig. 8c) changed from 278.7 K to 285.0 K by the single-stage approach, while from 278.4 K to 284.8 K for the two-stage approach.At M4 (Fig. 8d), the ground temperature by the single-stage approach and the two-stage approach changed from 274.4 K to 285.0 K and from 274.4 K to 284.8 K, respectively.
Owing to the small difference in ground temperatures at various locations calculated by the single-stage approach and the two-stage approach, no matter in 2D simulations or 3D simulations, it can be concluded that the two-stage approach can produce comparable results to the single-stage approach.Owing to the high computation efficiency, the two-stage approach is especially suitable for the simulation with a large domain size and long operation time.

Computation efficiency
All simulations in this study were carried out on a workstation with Windows 10 installed, an 8-core Intel(R) Xeon (R) W-2123 CPU running at 3.60 GHz, and 32.0 GB of RAM.Fig. 9 shows the computation time to complete 5-year-long 2D and 3D simulations by the single-stage approach and the two-stage approach.
As shown in Fig. 9a, under the identical computation settings as the single-stage approach, the two-stage approach took around 32% of the time to perform a 2D simulation.As shown in Fig. 9b, the two-stage approach took roughly 37% of the time as the single-stage approach for a 3D simulation with the same computation settings.
As mentioned earlier, the 2D and the 3D domains consist of 690 and 35,190 nodes, respectively.The computation time is proportional to the number of nodes in the domain.Based on Fig. 9, the computation time cost of a 2D simulation using the single-stage approach was about 1% of that of a 3D simulation.The computation time cost of a 2D simulation using the two-stage approach was around 0.9% of that of a 3D simulation.
Therefore, it can draw the conclusion that the two-stage approach outperforms the single-stage approach in terms of computation efficiency, and that 2D simulations are faster than 3D simulations.

Computation accuracy
Adopting the two-stage approach, the ground temperatures at the same buried depth obtained in 2D and 3D simulations are compared in Fig. 10.As shown in the figure, the ground temperatures from 2D and 3D simulations showed a similar development trend, while main difference can be found from January to April in a year.At the same buried depth, the ground temperatures obtained in the 2D simulation were lower than those in the 3D simulation during the first period of heating process.As the 2D simulation did not take the fluid temperature profile into consideration, a maximum difference ranged from 0.1 • C to 0.9 • C at various locations was caused.It should be mentioned that the time to complete a 2D simulation only accounted for 9% of the time to complete a 3D simulation under the same computation settings.
Although there are differences in the ground temperature obtained by the 2D simulation and the 3D simulation, 2D simulations by the twostage approach can be used for the preliminary design of the HGSHP system because of the high computational efficiency.

Conclusions
Numerical simulations to study ground thermal behaviour incorporating horizontal ground loops and realistic atmospheric processes are challenging and computationally expensive.In this paper, a newly proposed two-stage modelling approach is investigated for improving computational efficiency while maintaining accuracy of model predictions.
The model predictions of the two-stage approach were compared against the traditional single-stage modelling approach.A horizontal ground-loop with two 200-m-long legs, buried 3 m below the ground to extract heat, was evaluated using 2D and 3D simulations.The following conclusions can be drawn: • Good agreement is observed in the ground temperatures obtained by the two-stage approach and the single-stage approach.For instance, at the location of P3 in the 2D simulations, the ground temperature predicted by the single-stage approach and the two-stage approach ranged from 277.6 K to 284.5 K and from 277.6 K to 284.8 K, respectively.In 3D simulations, the ground temperature at the location of M3 calculated by the single-stage approach and the twostage approach varied from 278.7 K to 285.0 K and from 278.4 K to 284.8 K, respectively.• The two-stage approach is computationally robust.For example, in 2D simulations, the two-stage approach took around 32% of the time of the single-stage approach, while for 3D simulations it was roughly 37% of the time of the single-stage approach.• The two-stage approach can be used for modelling complex scenarios of realistic HGSHP systems installed in a large spatial domain and for long-term operation.For preliminary design of a HGSHP system, 2D simulations by the two-stage approach can be used considering the high computational efficiency and small differences in the ground temperature from 3D simulations.
The advantage of using the proposed two-staged approach would be more pronounced if the analyses are scaled up to more complex geometries and larger domains.In addition to computation time, computer energy consumption can be used as a metric to prove the benefit of developing highly efficient numerical methods and algorithms in the context of energy efficiency.
Finally, the numerical simulations in this paper were carried out based on available information from various sources, including monitored climatic condition, thermal loads, and soil parameters, further work can be done to validate the model results for an existing HGSHP system.

Declaration of Competing Interest
None.

Fig. 4 .
Fig. 4. Climatic variables: (a) ambient air temperature T air , (b) shortwave solar radiation H SW , (c) air relative humidity h air , (d) rainfall P, and (e) wind speed u y .

Table 2
Physical parameters of fluid in pipes.

Table 3
Ground surface boundary-related parameters.