Effect of groundwater forced seepage on heat transfer characteristics of borehole heat exchangers

A system of borehole heat exchangers (BHEs) combined with pumping–injection wells is established in areas where the groundwater is shallow and the seepage velocity is weak. The pumping and injection wells are set on both sides of the BHEs. According to the three-dimensional unsteady-state heat transfer model in the aquifer, the convection–dispersion analytical solution of excess temperature is derived that considers groundwater-forced seepage and thermal dispersion effects and axial effect of the BHEs. Then, we use the dimensional analysis method and similarity criteria to build a controllable forced seepage sandbox. The software FEFLOW 7.1 is adopted and the simulation results are validated by the theoretical analysis and the indoor experiment test. On this basis, the numerical simulation is used to explore the influence of different pumping–injection flow volume on the Darcy flow velocity of the aquifer where the BHEs are located, as well as the average heat transfer efficiency and the heat transfer rates with borehole depth. The results show that when the pumping flow volume increases from 200 m3 day−1 to 1200 m3 day−1, the Darcy velocity correspondingly increases to about 10 times. The average heat efficiency coefficient of the BHEs is increased by 11.5% in cooling stage, and by 7.5% in heating stage. When the pumping–injection flow volume is 400–600 m3 day−1, the increment of heat transfer rates of the BHEs reaches 12.8–17.9 W m−1 and 3.6–4.2 W m−1 per unit of borehole depth during the cooling stage and heating stage, respectively, and then decreases as the flow volume increases gradually.

part of the GHSP system, and the efficiency of the overall unit is directly affected by the magnitude of its heat exchange quantity. Therefore, enhancing the heat transfer performance of the BHEs is particularly important.
In general, the heat transfer performance of a BHE is studied using four methods: analytical method, numerical method, in situ thermal response test (TRT) and laboratory test. In most related fields, one or two of the above methods are adopted to study the effect of natural groundwater seepage on the heat transfer performance of the BHEs (Angelotti et al. 2014;Hu 2017;Choi et al. 2013;Wang et al. 2009).
Numerical simulation can achieve accurate solutions, which aid theoretical analysis. However, it needs extensive computational time. Hence, numerous calculations are often applied to the research and design of the BHEs, but the correctness and accuracy of the software should be verified beforehand (Casasso and Sethi 2014;Cuan et al. 2017). Hecht-Méndez et al. (2013) used FEFLOW, a finite element numerical simulation software, to perform the transient heat transport simulations for the fifteen scenarios of 25 boreholes and studied the distribution of the hydrodynamic field and the temperature field in a homogeneous confined aquifer.  analyzed the effect of unsaturated soil properties and groundwater flow on the performance of a GSHP system by the simulation software COMSOL Multiphysics. The numerical simulation was validated by combining experimental test results with the analytical solution that takes the multiplelayer substrates and groundwater flow into consideration.
Analytical solutions are preferred in most practical applications because of their excellent computational time and flexibility for parameterized design. The method of moving heat-source is applied in most analytical solutions to solve the problem of heat transfer under the influence of groundwater seepage (Li et al. 2015). Sutton and Nutter (2003) and Diao et al. (2004) presented an analytical solution considering groundwater flow. They both concluded that the temperature distribution in the vicinity of the boreholes was varied by the groundwater flow. Molina-Giraldo et al. (2011a) evaluated the influence of thermal dispersion on temperature plumes of geothermal systems by using analytical solutions. In the above analytical solutions, the axial effect is not considered and the borehole is considered to be a moving infinite line heat source. Molina-Giraldo et al. (2011b) proposed a moving finite line source model (MFLS) which considers the combined effect of groundwater flow and axial effects. However, the influence of thermal dispersion was not taken into account. Groundwater seepage affects heat transfer by involving gross heat convection and thermal dispersion, which is significant for the long-term temperature response of the BHEs (Luo et al. 2016). Therefore, an optimized analytical model, which considers the thermal dispersion effect in the MFLS model, is proposed in this paper.
TRTs applied to some real environments can truly reflect the heat transfer process of BHEs under the specific regional climatic characteristics and hydrogeological conditions. Most TRTs are performed on a test borehole to estimate the thermo-physical properties and borehole thermal resistance (Wang et al. 2009;Michopoulo et al. 2013;Choi and Ooka 2015). From the viewpoint of model validation, the data from TRTs are not suitable because TRTs' test times are short (usually 48-72 h) and susceptible to uncontrollable factors such as weather conditions during testing (Li et al. 2015;Smith and Elmore 2018;Shirazi and Bernier 2014). The indoor sandbox experiments enable the changing of some parameters, which is conducive to improve the speed and accuracy of the experiment (Zhang et al. 2015(Zhang et al. , 2016 because the experimental requirements have good controllability compared with the TRT. Most studies (Erol and François 2018;Samson et al. 2018;Hua et al. 2017;Zhang et al. 2014) point out that the presence of groundwater significantly affects the heat transfer between the BHE and its surrounding aquifer. However, most of these study conclusions are obtained under the condition of high natural seepage velocity, which is generally higher than 10 -8 m s −1 . Many areas such as the Bohai Rim plain have abundant groundwater reserves and the natural seepage velocity there is generally lower than 10 -8 m s −1 . In this case, it is usually considered that groundwater seepage has almost no effect on the heat transfer of the BHEs. Using the traditional GSHP system would render groundwater idle without having the positive effect that it should have. Therefore, it is necessary to optimize the traditional GSHP system under this special hydrogeological condition.
For the long-term application of GSHP systems, the cold and heat accumulation around the BHEs and the heat transfer quantity of the BHEs are decreased year by year. Scholars have proposed some optimized methods and improved measures to regulate the heat balance of the soil, such as solar-assisted GSHP systems (Østergaard et al. 2019;Nouri et al. 2019;Biglarian et al. 2019), and cooling-tower-assisted GSHP systems (Liu et al. 2018;Gong et al. 2018). Nouri et al. (2019) and Liu et al. (2018) study solar-assisted GSHP systems and cooling-tower-assisted GSHP systems, respectively. Their research results showed that the application of these hybrid systems could be helpful to reach considerable savings of energy through using free resources of stored heat in the ground, sun and air. However, solar-assisted GSHP systems and cooling-tower-assisted GSHP systems occupy a large land area and are greatly affected by seasonal and environmental factors.
Based on the local special hydrogeological conditions and existing GSHP system in Tianjin, China, this paper presents a BHE combined with pumping and injection wells located on both sides of the BHE. In reviewing the literature on the heat transfer performance of BHEs, we learn that most of them use one or two methods including analytical methods, numerical methods, in situ thermal response testing (TRT) and laboratory tests to study the effect of natural groundwater seepage on the heat transfer performance of the BHEs. So, we combine the theoretical analysis, experimental research, and the numerical simulation to discuss the changing trends of the Darcy flow velocity of aquifer, the dynamic changes of the average outlet water temperature, and average heat transfer efficiency of the BHEs. This paper aims to explore the influence of the different pumping-injection flow volume on the heat transfer characteristics of the BHEs under forced seepage. At the same time, it provides research ideas for reducing soil thermal accumulation, alleviating soil thermal and improving the heat transfer efficiency of the BHEs.

Geological survey
China's Bohai Rim region includes the plains of the four provinces of Liaoning, Hebei, Tianjin, Shandong and the surrounding mountainous areas (including the two peninsulas). The Bohai Rim plain was formed in the middle Pleistocene (Q2) and the late Pleistocene (Q3) to the Holocene (Q4) in the Quaternary period. The final formation of the Bohai Rim region is the result of the accumulation of sediment caused by rivers (Rong et al. 2012).
Tianjin is located in China's Bohai Rim plain, having abundant underground water and diverse hydrogeological structures. The whole Tianjin plain can be divided into the freshwater area, brackish water distribution underlying shallow freshwater area and brackish water area from north to south (Fig. 1). In the structure of aquifer, the sandy layer transforms into medium coarse sand, medium sand, medium fine sand, fine sand and silty sand from north to south.
The distribution area of shallow underground brackish water in Tianjin is 6922 km 2 , the brackish water area with a mineralization content of 2-3 g L −1 and 3-10 g L −1 are 3753 km 2 and 3169 km 2 , respectively, accounting for more than 2/3 of the city's total area (Zhou 2012). Groundwater resources are rich in reserves and convenient for exploitation, and the hydraulic gradient and the natural seepage generally range from 1.3 × 10 -2 m a −1 to 12 × 10 -1 m a −1 . Therefore, a BHEs combined with pumping-injection wells are suitable in the Tianjin Plain.

Underground thermo-physical properties
The GSHP system that was installed in Tianjin Binhai New Area, China, in 2016 is taken as the project prototype (Fig. 2). The project has a research area of 150 m × 120 m and a vertical depth of 83 m that it is divided into five geotechnical layers. Throughout the study area, a fine sand layer has good permeability and is regarded as a well-developed confined aquifer, which interacts with a silty sand layer to form a shallow-groundwater accumulation section. The geotechnical distribution and physical parameters are shown in Table 1.
The project, which is mainly responsible for the energy supply of the adjacent school, contains 219 BHEs with a center-to-center spacing of 4-5 m. The vertical depth of the borehole is 83 m and the 2U-DN32-HDPE-BHEs with a length of 80 m is arranged in the well. Expansive soil with sand is selected as the grout materials. Water is selected as the circulating refrigerant in the BHEs. The design parameters of the BHEs are shown in Table 2. Pumped well Inverted well Ground heat exchanger

Heat-seepage coupling model
It is assumed that the seepage process of the aquifer satisfies the following conditions: the physical parameters of aquifer and groundwater do not change with temperature; the seepage direction of groundwater is single and the vertical seepage process is ignored. According to the continuity equation of seepage flow and Darcy law (Nam et al. 2008), the continuity governing Eq. (1) and momentum governing Eq.
(2) in the anisotropic and homogeneous aquifers are established: In order to describe the change of heat transfer characteristics of the BHEs in the aquifer under complex space-time conditions accurately. The control equation of three-dimensional unsteady convection-thermal dispersion is established by taking the thermal dispersion effect into account, based on the convection-conduction equation (Capozza et al. 2013): in which the conductive part of the thermal dispersion tensor Λ cond and the dispersive part of thermodispersion tensor Λ disp are determined by Eqs. (4) and (5), respectively: The problem of determining the solution of seepage flow is associated with the problem of heat transfer by the momentum Eq. (2), so the heat-seepage coupling model is Volumetric heat capacity of refrigerant (water)/c r ρ r (J m −3 K −1 ) 4.18 × 10 6 Volumetric heat capacity of grout/c g ρ g (J (m −3 K −1 )) 2.19 × 10 6 Thermal conductivity of grout/λ g (W m −1 K −1 ) 1.90 constructed in the aquifer. The discharge (suction) heat process of the borehole can be considered as the source (sink) term in the aquifer thermal migration model. Because the diameter of the vertical borehole is much smaller than its depth, the heat transfer process of BHE is simplified to the heat transfer process of the moving finite line heat source in the semi-infinite medium.
On the basis of considering the thermal dispersion effect and the spatial position of the BHE, the moving finite line heat source model (MFLS) (Michopoulo et al. 2013) is improved. Under the premise of constant heat flow volume per unit length of BHE, the optimized analytical model (Eqs. 6, 7) is obtained by applying the method of images and the moving source theory. The optimized analytical model considers the heat convection, heat conduction and thermal dispersion which affects the process of groundwater seepage. In this way, the transient temperature in the aquifer caused by heat from the heat exchanger can be determined. Moreover, the model is coupled with the internal heat transfer process of the borehole.
In the analytic model, ρ e c e is the volumetric heat capacity of the porous medium (Eq. 8). When the thermal conductivity in the aquifers is the same in each direction, the thermal conductivity components of the Λ cond are determined by Eq. (9). λ x and λ y are the effective longitudinal and transverse thermal conductivity coefficients, which are determined by Eqs. (10) and (11). r is the distance to the source located on the z-axis at the (x 0 , y 0 , z') coordinates (Eq. 12): When the thermal dispersion effect in the aquifer is ignored (λ x = λ y= λ e ), the excess temperature analysis model (Eqs. 6, 7) can be simplified to Eqs. (13) and (14):

Experimental system
According to the engineering prototype, we design a complete laboratory experimental system to study a region, which is selected from a typical hydrogeological medium in a shallow aquifer, with the area of 20 m × 16 m and the buried depth of 24-44 m in Fig. 2. The engineering prototype contains six BHEs with the tube pitch of 4 m.

Design and construction of experimental system
The experimental system consists of a sandbox, heating apparatus, flow conditioning device and data-acquisition apparatus (Fig. 3). Since the seepage and heat transfer processes in the seepage sandbox and the engineering prototype follow the same form of Eqs. (1-3), the similitude relation ratio of basic design parameters is determined by dimensional analysis method according to the principle of similarity criteria (Li et al. 2015;Mao et al. 2006). The proportional relationship between hydrogeology and thermo-physical parameters of aquifers is 1:1 due to the use of equal volumetric weight filling of raw sand and seepage of raw water. In order to shorten the experimental period, the proportional relationship between the heat intensity of the BHEs and the difference of water head is set as 1:1. The geometric size proportional relationship between the actual and experimental was determined as 20:1. The Pr (Eq. 15), Fo (Eq. 16) and Pe (Eq. 17) in the engineering prototype and the sandbox are required to be equal in order to ensure that the experimental system can reproduce effectively the heat-seepage migration process of the aquifer. Thus, the operation time proportional relationship between the actual and the experimental is

Fig. 3
Schematic diagram of the experimental system determined to be 400:1 and the seepage velocity proportional relationship between the actual and the experimental is determined to be 1:20, so as to determine other design parameters, as shown in Table 3: According to the geometric similarity ratio, the sandbox is set to 1.2 m × 0.8 m × 1.2 m while the seepage region is 1 m × 0.8 m × 1 m and the confined aquifuge region is 1 m × 0.8 m × 0.2 m; the liquid supply/discharge region is 0.1 m × 0.8 m × 1.2 m, which is symmetrically set at both ends. Then, five overflow holes are drilled (diameter is 20 mm) with spacing of 0.2 m in the center line of the Plexiglas plate on the outsides of the liquid supply/discharge region. During the experiment, the hydraulic difference between the fluid supply and drainage areas is controlled by changing the opening and closing of the overflow holes at different heights, so as to change the seepage velocity in the seepage area. In addition, the holes (diameter is 5 mm) of the Plexiglas plate are evenly installed between liquid supply/discharge region and the seepage area, in order to ensure an even, horizontal seepage solution flow in the aquifer.
The K-type (± 0.1 ℃) thermocouple that has been treated with waterproof and anticorrosion treatment is selected to measure the temperature of the aquifer. Nine K-type thermocouples are embedded at 0.5 m from the bottom of the sandbox. Data acquisition equipment is used to record the temperature of the aquifer. Six 1.5-m-long PPR pipes, fixed on the bottom plate of the sandbox, are wound evenly with electric heating wire (50 W/m) to simulate the BHEs. The plane layout of the BHEs and temperature observation points are shown in Fig. 4.
The precision bath circulator THD-3015 is selected as the cold/heat source equipment to ensure that the temperature of the circulating groundwater meets the requirements of the test. Furthermore, the outside of the sandbox device, connecting pipes, storage tank and reservoir are pasted with thermal insulation materials with a thickness of 0.15 m. After that, six K-type thermocouples are fixed to all six sides of the sandbox to measure the heat loss of the experimental system.

Experimental scheme
The sandbox is filled with equal volumetric weight by layered wet filling method, and each layer is 50 mm thick and composed of raw sand. At the same time, the raw underground water is sprayed evenly and the sand layer is compacted to ensure that the porous medium in the sandbox has a unit weight of 1.68 ± 0.1 kg L −1 , which is similar to that of the underground aquifer. The raw sand is filled with a height of 1 m, and then a 0.2 m clay-gravel layer is laid on the upper part of the sandbox as a closed aquifer to isolate the aquifer from the external environment. The photograph of the sandbox without the thermal insulation material is shown in Fig. 5. After establishing the experimental system, central air-conditioning keeps the room temperature at 18 ℃. Raw water is injected below the surface at 16 ℃ into the sandbox. When all temperature measuring points in the sandbox are maintained at 16 ± 0.1 ℃, the water level in the piezometer is stable and no bubbles appeared, the porous medium of the sandbox can be considered to be saturated. Before the experiment began, it is proved that Darcy velocity was 2.4 ± 0.02 × 10 -6 m s −1 (Pe ≈ 2.5) when the hydraulic difference is 0.2 m, which satisfies the proportional relationship between actual and experimental seepage velocity. During the experiment, the 16 ℃ underground raw water with 0.2 m hydraulic difference is filled continuously into the sandbox to ensure a stable Darcy velocity. Simultaneously, the BHEs 1-4 are opened in the test. The experimental running time is set at 5.4 h and the temperature of each measuring point is recorded at 1-min intervals. Through the experiment, it is found that the temperature measurement range of the six thermocouples installed outside the sandbox is 18 ± 0.5 ℃. Therefore, the heat loss from the sandbox is negligible.

Compared with analytical and experimental data of the excess temperature
According to the engineering prototype, FEFLOW 7.1 is applied to establish the geometric model, conduct the mesh division (triangular non-isometric) and set the parameters, and then the numerical simulation is performed. Meanwhile, MATLAB 2012 is used to calculate the transient temperature response caused by the running BHE, according to the unsteady analytical model of the excess temperature in the aquifer (Eqs. 6, 7). Then, the experimental results and the analytical solutions are compared with the numerical solution.
In the simulation, the clay layer and the fine sand layer are confined aquifers, while the thickness of fine sand layer is 20 m. The clay layer in the upper part of the study area is defined as an impervious and adiabatic boundary, ignoring the influence of external environment. The four flanks of the 20 × 16 × 83 m 3 model are defined as the fixed hydraulic head and constant temperature boundary. Then, the parameters of the six BHEs in the region are set according to Table 2. In addition, the triangular non-isometric mesh is adopted in the entire geotechnical layer. The total number of nodes in the physical model is 51740, and the number of grids is 93156. The fixed-time step method is used to solve the problem. When the time-step length is set to 1 min, there is a total of 324 time steps, and the maximum iteration is 2500 times per step.
For facilitating the comparison, the time t, the excess temperature ΔT and the coordinate displacement x are all in a dimensionless form. Due to the different proportional relationships between the engineering prototype and the experimental system in operation time and geometric size, Fo is taken as dimensionless time (Eq. 16), θ is taken as dimensionless excess temperature (Eq. 18) and L is taken as dimensionless coordinate displacement (Eq. 19). Meanwhile, the root mean square error (RMSE) of dimensionless excess temperature θ (Eq. 20) is selected as the similarity index between experimental results, analytical solutions and numerical solutions: in which θ N(i)/A(i) θ E(i) correspond to the numerical solutions, analytical solutions and experimental result, respectively. The RMSE between the analytical solutions and the experimental result in Fig. 6a is < 5%, while the RMSE between the numerical solution and the experimental result is < 1%. Besides, the RMSE between the analytical solution and the experimental result is 3.8%, and the RMSE of the numerical solutions and experimental results is 0.5% in Fig. 6b. The results show that the analytical solution and experimental results are consistent with the numerical solution. So, FEFLOW 7.1 can simulate the heat transfer process of the BHEs in the aquifer.

Compared with existing studied data of the excess temperature
The numerical solutions regarding the mean temperature change of the aquifer are then compared to available data derived from Hecht-Méndez et al. (2013). Since the studies refer to different ground properties, groundwater velocities, and dimensions, the comparison is performed by calculating for each condition in the papers the Pe, Fo, θ (Table 4).
It is found that the "scenario 12" in Hecht-Méndez's research (Hecht-Méndez et al. 2013) is close to the working conditions in this paper through the comparative research of Pe (Fig. 7). The results show that the maximum error is 10.3% and other errors are < 10%. Therefore, FEFLOW 7.1 is effective and correct to simulate the heat transfer process of the BHEs under the forced seepage.

Analysis of examples
This section mainly explores the influence of the pumping-injection flow volume of the BHEs combined with pumping-injection well on enhancing the heat transfer effect, compared with the established traditional GSHP system.

Generalized model and operation scheme
Due to the complex and diverse layout of the on-site well group and large number   of boreholes, the BHEs (Fig. 2) are generalized into five groups of the BHEs in Fig. 8. There are seven boreholes in each group, and the spacing of each borehole is 4-5 m. To ensure the groundwater synchronous recharge, two pumping wells and five injection wells are set. The design parameters of BHE are shown in Table 2. The seven working wells are all incomplete diving wells with a depth of 60 m and a diameter of 0.4 m. Meanwhile, the filter is placed at a depth of 30-40 m.
Based on the well group layout, the horizontal calculation area is 150 m × 120 m and the vertical calculation range is 0-83 m. From top to bottom, the rock-soil layer with a thickness of 83 m is separated into five horizontal fault types. The spatial distribution and physical properties of this layer are shown in Table 1.
In the numerical simulation, the effects of precipitation and evaporation are ignored, the silty clay layer in the upper and lower parts of the study area is defined as the impervious and adiabatic boundary. The four flanks of the physical model are defined as the fixed hydraulic head and constant temperature boundary. In addition, the pumping-injection wells are defined as the constant flow boundary. Each soil layer is separated by a triangular non-isometric mesh. Local refinements of the mesh are set at the BHEs and the pumping-injection well locations.
The total number of nodes and grids in the physical model is 248,752 and 413,295, respectively. The fixed-time step method is used to solve the problem. The time step is set to 1d, the total time step is set to 3650, and the maximum number of iterations is set to 2500 (times per step). According to the established model, the simulation calculation of eight operating modes is carried out (Table 5).
One operation cycle (1 year) is divided into 4 stages, which are summer cooling stage (4 months), autumn intermittent stage 1 (2 months), winter heating stage (4 months) and spring intermittent stage 2 (2 months). The system runs five operation cycles and the BHEs running continuously for 10 h per day during the cooling and heating stages. The BHEs are set to unconstant power heat absorption (dissipation) operation. The inlet water temperature of the BHEs during cooling/heating stage is constant at 31 ℃/6 ℃, respectively. The temperature of the outlet water of the BHEs and the temperature of the inlet and outlet water on the side of the heat pump unit change gradually with the operation of the system.

Analysis and discussion
The Darcy velocity and hydrodynamic distribution are significantly different in the same aquifer due to the difference in the total pumping-injection flow volume (∑G) in Fig. 9. Darcy velocity, which increases approximately tenfold, increases from 2.4-3.2 × 10 -7 m s −1 to 2.0-3.0 × 10 -6 m s −1 when the total pumping-injection flow volume increases from 200 m 3 d −1 to 1200 m 3 d −1 . Due to difference pumping-injection flow volume, the temperature fields of the same aquifer profile have a significant difference (Fig. 10). Taking well group E of the BHEs as the research object, the injection well and the pumping well with an interval of 120 m are defined as the upstream boundary and downstream boundary, respectively. To describe accurately the evolution process of the aquifer's temperature field under different operation modes, the calculation area with a temperature change of ± 0.5 ℃ is defined as the thermal diffusion range of the BHEs. Besides, the heatinfluencing radius is defined as the coordinate distance between E (g) and the farthest acting position.
When the pumping-injection flow volume is 0 m 3 d −1 , there is only the heat conduction between the BHEs, aquifer and aqueous medium units. The heat transfer process is slow and the heat influence range is diffused symmetrically around the BHEs. By the end of the cooling stage (120 d), the thermal radius of the BHEs in both upstream and downstream areas is 11 m.
The range of thermal diffusion in the downstream area of the BHEs significantly expands with the pumping-injection flow volume increases from 200 m 3 d −1 to 1200 m 3 d −1 . At 120 d, the thermal radius along the direction of pumping reaches 19 m, 35 m and 49 m, so the thermal radius of 1200 m 3 d −1 is 2.6 times of that of 200 m 3 d −1 . With the increase of the flow volume, the migration speed of the temperature fronts accelerates and the thermal radius enlarges continuously; moreover, the thermal radius in the upstream zone is smaller than 11 m.
According to the theory of heat and mass transfer, it is precisely because of the increase of Darcy velocity that the convective heat transfer intensity and thermomechanical dispersion effect are improved correspondingly, thereby expanding the range of thermal diffusion in the downstream region and alleviating the thermal accumulation phenomenon of BHE. So, the dynamic variation of ∆T disp /∆T cond with time is calculated (Fig. 11) in order to obtain the difference between the temperature response ∆T disp with forced groundwater seepage and the temperature response ∆T cond without groundwater seepage at different temperature observation points.
As shown in Fig. 11a, the temperature is decreases gradually. When the running time exceeds 50% of the whole period (Fo ≥ 0. 2), ∆T disp /∆T cond < 1. As shown in Fig. 11b, ∆T disp /∆T cond tends to be stable. At the end of the cooling stage (120 d), ∆T disp /∆T cond > 1. The change of Darcy velocity has no obvious influence on the change rate and amplitude of excess temperature because the thermal convection and thermal dispersion enhances thermal interference between the BHEs. As shown in Fig. 11c, the excess temperature tends to be stable and the curve of ∆T disp /∆T cond tends to be smooth. The mutative degree of excess temperature enhances with the increase of Darcy velocity. However, the difference between ∆T cond and ∆T disp tends to decrease when Darcy velocity increases to a certain extent.
In this paper, the ratio of actual average heat exchange quantity ( Q ) to theoretical average heat exchange quantity ( Q ) of the BHEs is defined as the average heat transfer efficiency coefficient ( Ē ). The degree to which the BHEs performance deviates from the theoretical heat transfer under different operating conditions is revealed. At the early stage of various operation modes, due to the large temperature difference in heat transfer between the boreholes and the rock-soil layer, the temperature difference between the inlet and outlet water of the BHEs is higher and the corresponding energy efficiency coefficient also rises (Fig. 12). The heat exchange rate is reduced because of the decreasing temperature difference between the boreholes and the surrounding medium. As a result, the temperature difference of the inlet and outlet water of the BHEs decreases rapidly. At the end of the fifth operation cycle, the average heat transfer efficiency of the five types of operation modes in the cooling stage is 7.1%, 8.5%, 11.1%, 13.4% and 18.6% while in the heating stage is 10.6%, 11.9%, 12.8%, 14.6% and 18.1%, respectively.
The cumulative distribution curve can be used as another evaluation index to describe the duration of a certain heat exchange efficiency of the BHEs (Fig. 13). The cumulative distribution curve means that under cooling and heating stages, the absolute value of the water temperature difference between the inlet and outlet of the BHEs corresponding to the wellbore during the operation time is counted and arranged in descending order, and the temperature difference between the inlet and outlet of the BHEs is a cumulatively time-varying distribution curve.
The absolute values of the inlet and outlet temperature difference of the BHEs is arranged in descending order, and the cumulative average temperature difference distribution in the cooling stages and heating stages over the whole simulation numerical stages (5 years) is calculated. When the cumulative time exceeds 20% of the total operation stage, the temperature difference tends to be moderate. When the cumulative time reaches 50%, and the temperature difference reaches a steady state.
The inlet and outlet temperature difference of the BHEs corresponding to the median time is used as another evaluation index to compare the heat transfer performance of the BHEs under different modes. In the cooling stage, the median of temperature difference is 2.34 ℃ when the pumping-injection flow volume is 0 m 3 day −1 . When the pumping-injection flow volume rises from 200 m 3 day −1 to 1200 m 3 day −1 , the growth rate of the median temperature difference varies from 11.5% to 73.9%. In the heating stage, the median temperature difference is 1.76 ℃ when the BHEs run individually. When the pumping-injection flow volume increases from 200 m 3 d −1 to 1200 m 3 d −1 , the growth rate of the median temperature difference increases from 10.2% to 34.1%. With the increase of pumping-injection flow volume, the Darcy velocity and the heat exchange intensity, the decline rate of the temperature difference slows down while the time required reaching a steady-state increases.
The heat transfer per unit depth of the BHEs ( q ) is determined by the circulating water volume, the thickness of different rock and soil layers, and the temperature difference between the inlet and outlet water. Therefore, under different operating conditions, q is used as an evaluation parameter for the heat transfer characteristics of the BHEs in various rock and soil layers. The q rises gradually when pumping and injection flow volume increases (Fig. 14). Taking the fifth year as an example, the q in cooling and heating stages with the increase of the pumping-injection flow volume from 0 m 3 d −1 to 1200 m 3 d −1 varies from 30.3 W m −1 and 25.4 W m −1 to 82.2 W m −1 and 39.6 W m −1 , respectively.
Thermal accumulation becomes serious with the continuous operation of the system, the average heat transfer efficiency gets lower on the condition that the pumping Fig. 13 Cumulative distributions of the temperature difference of the inlet and outlet water of the BHEs: a the cooling stage; b the heating stage and injection wells cease to work or the pumping-injection flow volume becomes smaller (200 m 3 d −1 ). As a result, cooling and heating quantity declines every year. However, the rock-soil layer is not a single cold or heat source. An energy storage body stores and emits heat to a certain extent. Therefore, the heat storage capacity of the rock-soil layer is enhanced by thermal accumulation in the wells group's area across the season in the summer. During the heating stage, the temperature of the rock-soil layer is higher than the initial stage. The temperature difference between the soil and the circulating solution increases so that the heat exchange quantity of the BHEs is improved at this stage.
Under the cooling and heating stages, the relation between Δ q (known as the increment of average heat transfer rate per unit depth of BHEs) and Pe is determined by fitting formula which shows that Δ q is distributed as a Gaussian function with the uptrend of Pe (Fig. 15). Although the effect of convective heat transfer and thermal dispersion between the BHEs and the rock-soil layer can be enhanced effectively by strengthening the velocity of groundwater seepage, the average heat transfer rates per unit borehole depth of the BHEs (q) does not raise linearly with the increase of the pumping-injection flow volume.
The research shows that Darcy velocity is only 0.6 × 10 -6 -1.4 × 10 -6 m s −1 when the pumping and injection flow volume is 400-600 m 3 d −1 , but the Δ q reaches 12.8-17.9 W m −1 and 3.6-4.2 W m −1 during the cooling stage and heating stage which are located on both sides of the extremum point of the distribution curve. As the designed pumping-injection flow volume further increases, not only the Δ q decreases gradually, but also the energy consumption of pumping and injection pumps increase which leads to the increase of operational cost of the system. Furthermore, the change of aquifer spatial structure will be irreversible if the forced seepage velocity is too high. Therefore, in order to obtain the best heat transfer enhancement effect, system environment and economic benefits, the best reference range for the BHEs is taken to be when Darcy velocity reaches 0.6 × 10 -6 -1.4 × 10 -6 m s −1 . So, the best flow volume range of pumping-injection wells is 400-600 m 3 day −1 in this paper.

Conclusion
In this article, a system of BHEs combined with pumping-injection wells is presented. Based upon the three-dimensional unsteady-state heat transfer model in aquifer, a convection-dispersion analytical solution of excess temperature in aquifer is derived that considered groundwater-forced seepage and thermal dispersion effects in aquifer and the axial effect of the BHEs. At the same time, according to the dimensional analysis method and the similarity criteria, a forced seepage indoor 3-D sandbox is established. Then, FEFLOW 7.1 is used to discuss how different pumping-injection flow volumes affect the heat transfer characteristics of the BHEs. The main conclusions of this study can be summarized as follows: Fig. 15 The total energy consumption of pumps and the increment of heat transfer rate per unit depth of the BHEs (Δq = q t+1 − q t ) change over Pe (a) Based on the established mathematical model, an indoor sandbox is built and an analytical solution is developed. According to the engineering prototype, FEFLOW 7.1 is adopted to establish a geometric model for the numerical simulation. The results indicate that the analytical solutions and experimental results are roughly consistent with the numerical solutions. The RMSE between the analytical result and the experimental result is 3.8%; the RMSE between the numerical result and experimental results is 0.5%. Thus, the numerical simulation software FEFLOW 7.1 can simulate effectively and correctly the evolution process of the aquifer temperature field during the taking heat of the BHEs. (b) When the pumping-injection flow volume increases from 200 m 3 day −1 to 1200 m 3 day −1 , Darcy velocity in the fine sand layer increases from 2.4-3.2 × 10 -7 m s −1 to 2-3 × 10 -6 m s −1 and the effect of convective heat transfer and thermal dispersion is intensified. In the downstream region of the BHEs, the thermal action radius is 2.6 times of the original radius. Hence, the BHEs combined with pumping-injection well can reduce soil thermal accumulation and alleviate soil thermal interference effectively. (c) For long-term running of the BHEs, the cumulative distribution curve of the inlet and outlet water average temperature difference is introduced to describe the duration of a certain heat exchange efficiency of the BHEs. The average temperature difference corresponding to the median time is used as a new parameter for evaluating the heat transfer performance of BHEs. In the cooling stage, the pumping-injection flow volume rises from 200 m 3 d −1 to 1200 m 3 d −1 , the growth rate of the median temperature difference varies 62.4%. In the heating stage, the pumping-injection flow volume increases from 200 m 3 d −1 to 1200 m 3 d −1 , the growth rate of the median temperature difference increases 23.9%. At the same time, the average energy efficiency coefficient of the BHEs improves 11.5% (cooling stage) and 7.5% (heating stage). (d) The simulation indicates that in cooling and heating stages, the heat transfer characteristic of the BHEs increases significantly with the increase of the Darcy velocity. However, the relation curve between the Δ‾q and Pe has a Gaussian function distribution in cooling and heating stages. Therefore, the heat exchange quantity of the BHEs cannot be increased continuously by increasing the pumping-injection flow volume infinitely. The pumping-injection flow volume that can make Darcy velocity reach 0.6 × 10 -6 -1.4 × 10 -6 m s −1 is the best reference range of the BHEs combined with pumping-injection well. (e) The default geographic parameters are the homogeneous of the same aquifer, and the thickness is 20 m. However, the greater the thickness of the aquifer, the better the heat transfer characteristics of the BHEs. The influence of the thickness on the heat transfer characteristics is not the focus of this article. Later, the influence of aquifer thermo-physical properties (thermal conductivity, porosity, permeability, etc.) and heterogeneity on the heat transfer characteristics of the BHEs will be studied. Greek symbols μ: Dynamic viscosity (Pa s); ε: Porosity; Λ: Tensor of thermal hydrodynamic dispersion (W m −1 K −1 ); δ ij : Kronecker tensor (i = j, δ ij = 1; i ≠ j, δ ij = 0); γ: Compressibility coefficient; λ: Thermal conductivity (W m −1 K −1 ); α L , α T : Longitudinal and transverse thermal dispersivity, respectively, of fluid (m); ρc: Volumetric heat capacity (J m −3 K −1 ); θ: Dimensionless temperature.