Factors Affecting the Installation Potential of Ground Source Heat Pump Systems : A Comparative Study for the Sendai Plain and Aizu Basin , Japan

Evaluating the installation potential of ground source heat pump (GSHP) systems based on the hydrogeological condition of an area is important for the installation and sustainable use of the system. This work is the first to have compared the distributions of heat exchange rate in the Sendai Plain and Aizu Basin (Japan) in terms of topographical and hydrogeological conditions. A regional groundwater flow and heat transport model was constructed for the Sendai Plain. Suitability assessment was conducted for an identical closed-loop system by preparing the distribution maps of heat exchange rate for space heating for the plain and basin. For both locations, the upstream area showed a higher heat exchange rate than the downstream area. Multiple regression analysis was conducted using heat exchange rate as a response variable. Average groundwater flow velocity and average subsurface temperature were considered as explanatory variables. The heat exchange rate for the plain, whose Péclet number ranged from 3.5 × 10−3–7.3 × 10−2, was affected by groundwater flow velocity and subsurface temperature. The exchange rate for the basin, whose Péclet number ranged from 8.5 × 10−2–5.8 × 10−1, was affected by groundwater flow velocity. Inland basins are likely to be more suitable for GSHP system installation utilizing groundwater flow than coastal plains in terms of inclination of slope. This study showed that multiple regression analysis can reveal factors affecting the heat exchange rate as well as the degree to which they affect it.


Introduction
Ground source heat pump (GSHP) systems are an energy efficient and environment friendly technology.They use natural heat energy of the subsurface stored at the shallow depths for space heating and space cooling among other uses [1,2].GSHP systems have been installed since 1980s and numerous evaluations of these systems have been performed [3][4][5].The development and uptake of this system in Japan are gradually increasing [6].However, the number of installations of GSHP systems is lower than in other locations around the world such as America and Europe.This is mainly because it is expensive to install such systems, particularly when the installation capacity exceeds the required capacity.Thus, it is necessary to identify suitable areas to install GSHP systems, to ensure their appropriate design and sustainable growth.In Japan, GSHP suitability mapping has been carried out for the Tsukushi Plain [1] and the Tsugaru Plain [7] based on hydrological parameters such as groundwater flow, groundwater temperature, and hydraulic head, because GSHP efficiency is affected by groundwater flow.Subsurface temperature is affected not only by heat conduction but also by heat advection of groundwater flow [6][7][8][9][10][11]. Therefore, the groundwater condition of the region influences the heat exchange rate and performance of GSHP systems [12][13][14][15].Major cities and towns in Japan are located on plains and basins that are part of the Quaternary system, where groundwater flows actively [16].Therefore, it is important to design the GSHP system considering groundwater flow.
Various methods have been used to evaluate installation potential of the systems in Japan.For example, Uchida et al. [17] prepared four kinds of thematic maps using the results of a numerical model and hydrogeological information (groundwater flow velocity, subsurface temperature, depth of groundwater table, and ratio of sand-gravel layers).Ohtani et al. [18] examined the distribution of groundwater flow velocity based on geological and hydrogeological information using geographical information systems to evaluate the installation potential of a GSHP system.Shrestha et al. [19] evaluated the suitability of a GSHP system in the Aizu Basin by creating distribution maps of heat exchange rate.Topographical and hydrogeological conditions differ between plains and basins, as does groundwater flow, which affects the heat exchange rate.However, no study has compared the GSHP system installation potential of one region with that of another in terms of these conditions.It is not clear whether and how differences in these conditions between basins and plains may affect the heat exchange rate.Comparing these regions under the same computational conditions for the system will shed light on the factors affecting the heat exchange rate, and thereby, the installation potential of the system.
The main objectives of this study are to compare the distributions of heat exchange rate in the Sendai Plain and Aizu Basin in Japan and to examine factors affecting this distribution owing to differences in topographical and hydrogeological conditions.Thus, this study calculates the heat exchange rate in the Sendai Plain, assuming the same computational conditions as those used for the Aizu Basin in a previous work [19] and analyzes and compares the situation of the plain with that of the basin.In this study, 3D groundwater flow and heat transport model were constructed for the Sendai Plain, and the model was verified based on field survey data.Next, a suitability assessment was conducted for a general closed-loop system by preparing a distribution map of heat exchange rates for space heating.Following Shrestha et al. [19], this map was prepared by using the results of the heat exchange simulations conducted with ground heat exchanger (GHE) models constructed at several locations within the Sendai Plain.Furthermore, the heat exchange rate distribution map for the Sendai Plain was compared with that of the Aizu Basin [19].This study contributes to the literature by comparing the distributions of heat exchange rate in the two regions in terms of topographical and hydrogeological conditions and by examining factors affecting the rate based on multiple regression analysis.Given the need for GSHP system, this approach is valuable for examining the factors affecting the heat exchange rate in various plains and basins in Japan and other countries.

Study Area
The Sendai Plain, which is located in Miyagi Prefecture, stretches in the N-S and W-E direction over a distance of 40 km and 10 km, respectively (Figure 1).The maximum thickness of the Quaternary formation in the plain is about 80 m [20].Small-scale alluvial fans exist in the inland area, and low-lying shore and dunes are observed in the seacoast area in the plain.The upper boundary depth of the gravel layer tends to increase from the inland area to the seacoast area.Mud and sand layers are underlain by gravel around the seacoast area.

3D Groundwater Flow and Heat Transport Model (Regional Model)
The 3D groundwater flow and heat transport model were constructed using finite element software FEFLOW [21] to evaluate the suitability of the GSHP system.The governing equation for groundwater flow and heat transport simulation in the saturated zone with variable-density fluid is given by: where S 0 is the specific storage [L −1 ], h is the hydraulic head [L], t is the time [T], q is the Darcy velocity The size of the modeled area was about 3600 km 2 , considering the surface-water divide in the Sendai Plain (Figure 2).The horizontal mesh of the model was refined along the rivers and major areas of the plain.In the vertical mesh, layers 1-2 belong to the Quaternary system, layers 3-9 belong to the upper Neogene, and layers 10-17 belong to the lower Neogene.The Quaternary system and upper Neogene act as the main aquifers of the plain, where groundwater flow primarily occurs.Lower Neogene was regarded as the hydrogeological basement of the Sendai Plain.The thickness of the Quaternary system, upper Neogene, and lower Neogene was up to 80 m, 540 to 620 m, and 300 m, respectively.Table 1 shows the hydrogeological and physical parameters of the model.These parameters were based on Bair and Lahm [22] and the Geo-Heat Promotion Association of Japan [23].Volumetric heat capacity of 2.52 MJ/m 3 •K was assigned to the entire model.Several boundary conditions were set to conduct groundwater flow and heat transport simulation with the model.Constant head boundaries were assigned at the model's top boundary to analyze groundwater flow.The hydraulic heads at the constant head boundaries outside the plain were assigned using a river level map prepared for a 1 km × 1 km mesh using digital elevation model data sourced from the Geospatial Information Authority of Japan (GSI).The model's bottom and vertical sides were regarded as impermeable layers.Constant temperature boundaries were set at the model's top surface and bottom surface to analyze the subsurface thermal regime.Average annual air temperature of the location was assigned as a prescribed value at the constant-temperature boundary at the model's top surface, considering the variation in elevation.On the other hand, the constant temperature value at the model's bottom surface was assigned after estimating the temperature calculated using the vertical temperature gradient of 0.035 • C/m, which is the typical value for Japan.Groundwater in the Sendai Plain is mainly used for industrial and agricultural purposes, and the corresponding values were estimated and inputted in the model as groundwater discharge points after considering unpublished survey reports from Miyagi Prefecture.

Validation of the Constructed Model
The constructed model was validated by comparing the calculated values of hydraulic heads, subsurface temperature profiles (vertical distribution of subsurface temperature), and thermal response test (TRT) with their corresponding measured field data for the groundwater monitoring wells in the plain.Figure 3 shows location of the groundwater monitoring wells.Figures 4-6 show the comparison of measured hydraulic heads, subsurface temperature profiles, and TRT with their calculated data.The data on hydraulic heads were sourced from unpublished survey reports of Miyagi Prefecture, which measured nine wells located in the plain.Subsurface temperature profiles were reported by Uchida et al. [24], who measured the profiles in the plain.TRT data were sourced from Kaneko et al. [25], who surveyed urban areas of the plain.The calculated values were higher than the measured values for some wells and the maximum difference between the two was about 3 m at site S-9, which contains four groundwater observation wells (Figure 4).This is because an industrial complex exists around the site S-9 and groundwater may be pumped out to meet the complex's water requirement.However, the locations of the pumping wells and their screen depths were not known.More accurate estimation of such information, such as location, pumping rates, and screen depth, can improve this model.The coefficient of determination between the calculated hydraulic heads and measured values was 0.71 in this analysis and the tendency of calculated hydraulic head was reproduced.A comparison of the calculated and measured vertical subsurface temperature profiles is shown in Figure 5.The calculated values for sites W-1 and R-1 did not agree with the measured ones.These sites are located in the groundwater recharge area of the plain.Some faults exist in the plain.They influence the accuracy of the estimated thickness of the Quaternary layer.The calculated values for the other wells showed good agreement with the measured values or reproduced temperature gradient of the profile.TRT was carried out at a Borehole Heat Exchanger (BHE; well depth: 101.35 m) in Sendai's urban area from 4th to 11th September 2017 (Figure 3).The heated water was circulated through the U-tube for about 66 h.The heat exchange rate between the formation and heat exchanger was 4.687 kW.Inlet temperature was used as an input parameter and outlet temperature was used as the matching parameter of the constructed model.With regard to the verification of the TRT, the calculated values agreed with the measured values, and the mean absolute error between the two was 0.33 (Figure 6).Therefore, it is logical to say that the calculated groundwater flow and subsurface thermal regime of the constructed model are consistent with the actual groundwater conditions in the Sendai Plain.

Ground Heat Exchanger Model
A suitability map showing the distribution of heat exchange rates is useful for examining suitable areas for the installation of a GSHP system.Assuming a general closed-loop system, identical ground heat exchanger (GHE) models of dimensions 20 m × 20 m × 120 m were constructed at 33 locations (Figure 7) to calculate heat exchange rates at these locations.These analyses were conducted for space heating purposes, to evaluate the heat exchange rate under the same conditions as those in Shrestha et al. [19].
GHE models were constructed using FEFLOW [20].Parameters inputted in the GHE models were the same as those assigned to the regional model, as shown in Table 1.The same geological data and hydrological and thermal parameters from the same GHE locations in the regional model were used for each GHE model.Similarly, initial and boundary conditions of the GHE models were inputted based on the hydraulic head, groundwater flow velocity, and subsurface temperatures calculated from the regional model at each GHE location.
A GHE was set at 100 m depth at the center nodes of the GHE models.A double U-tube of diameter 34 mm and grout of silica sand were considered.
Several boundary conditions were set to conduct groundwater flow and heat transport simulations for the GHE models.The models' top and bottom were regarded as impermeable layers.Constant head boundaries were assigned at vertical faces of the GHE models to obtain the same flow velocities, which were calculated at the same points as those in the regional model.Constant temperature boundaries were assigned to all layers based on results obtained from the regional model.The settings for the models were the same as those in Shrestha et al. [19].

Distribution Map of Heat Exchange Rates
Heat exchange simulations were conducted for space heating for each GHE model at the respective locations.The operating scenario was set as 120 d of space heating per year from December to March assuming 24-h operation.The inlet temperature and flow rate of the circulation fluid were set as 5 • C and 20 L/min, respectively.Water was used as the circulation fluid.This setting of GHE operation was the same as that in Shrestha et al. [19].The heat exchange rates at the GHE locations were estimated from the simulations.Figure 8 shows the distribution map of heat exchange rates for the Sendai Plain.The average, maximum, and minimum values of heat exchange rates were 22.5, 38.6, and 16.3 W/m, respectively.Heat exchange rates were higher at the periphery of the plain and lower in its coastal areas.Figure 9 shows the distribution map of average groundwater flow velocity from the surface to 100 m depth.These values were based on the results of the regional model.Groundwater flow velocity was also greater along the periphery of the plain and lower in its coastal areas.Greater groundwater flow velocity increases the apparent thermal conductivity of the subsurface, which is caused by heat transfer through advection of groundwater flow [7,[26][27][28].Therefore, the greater velocity may also help retain the thermal energy at the GHE by keeping a constant subsurface temperature and improving the heat exchange rate of the system at and around the GHE location.
The highest heat exchange rate (38.6 W/m) was observed in the terrace area of the Sendai Plain.This site may be largely affected by subsurface temperature rather than groundwater flow.Figure 10 shows the distribution map of the average subsurface temperature from the surface to 100 m depth.The site with the highest heat exchange rate showed the highest average temperature of the 33 GHE locations.Space heating operation causes a decrease in subsurface temperature.Higher average subsurface temperature may increase the heat exchange rate, because the temperature difference between the inlet temperature of the GSHP system and the subsurface temperature increases.

Comparison of Heat Exchange Rates in the Sendai Plain with in the Aizu Basin
Figure 11 shows the distribution map of the heat exchange rates in the Aizu Basin.The average, maximum, and minimum values of heat exchange rates were 35.6, 44.7, and 28.6 W/m, respectively.Groundwater flow velocity was higher in the northern and southern areas of the basin.This is because these areas, which are located in the upstream area of the basin, have higher hydraulic gradients and geology amenable to higher velocities [19].This tendency of the heat exchange rate being higher in the upstream area compared to its downstream counterpart was also observed in the Sendai Plain.
Comparing the heat exchange rates in the Sendai Plain and Aizu Basin shows that the heat exchange rates for the former are lower.This is because groundwater flow velocity through the borehole heat exchangers (BHEs) in the Sendai Plain is lower than that in the Aizu Basin. Figure 12 shows the results of the regression analysis.Based on the simulation results of the GHE models for both the Sendai Plain and Aizu Basin, the regression analysis was carried out using the heat exchange rate as a response variable and average groundwater flow velocity or average subsurface temperature as the explanatory variable (Figure 12a-d).Average groundwater flow velocity and average subsurface temperature in the Sendai Plain are depicted in Figures 9 and 10, respectively.Average groundwater flow velocity and average subsurface temperature in the Aizu Basin were based on the simulation results reported by Shrestha et al. [19].The Aizu Basin has about 10 times higher groundwater flow velocity than the Sendai Plain.The coefficients of determination between the average groundwater flow velocity and heat exchange rate for the plain and basin were 0.27 and 0.90, respectively (Figure 12a,b).The average subsurface temperature of the plain ranged from 12 to 21 • C, and the corresponding value for the basin ranged from 12 to 16 • C. The coefficients of determination between the average subsurface temperature and heat exchange rate for the plain and basin were 0.24 and 0.22, respectively (Figure 12c,d).Multiple regression analysis was conducted using the following equation: where Y is the heat exchange rate (response variable); X 1 and X 2 denote the average groundwater flow velocity and average subsurface temperature at a depth of 0-100 m, respectively (explanatory variable); and a, b, and c are coefficients.The results of the regression analysis are shown in Table 2.The adjusted R-square for the plain and basin was 0.89 and 0.77, respectively (Figure 12e,f).For the plain, the t test for both X 1 (average groundwater flow velocity) and X 2 (average subsurface temperature) provided a result greater than 2, and the p values of both X 1 and X 2 were smaller than 0.05.For the basin, the t test for X 1 provided a result greater than 2, and the p value of X 1 was smaller than 0.05.However, the result of the t test for X 2 was smaller than 2, and the corresponding p value was greater than 0.05.These results show that the heat exchange rate was related to groundwater flow velocity and subsurface temperature for the plain, and only groundwater flow velocity for the basin.

Discussion
The heat exchange rate was related to groundwater flow velocity and subsurface temperature for the Sendai Plain, and only groundwater flow velocity for the Aizu Basin.This difference may be attributed to groundwater flow velocity.Whereas average groundwater flow velocity at the BHEs in the Sendai Plain range from 6.3 × 10 −4 to 1.3 × 10 −2 m/d, the groundwater flow velocity of the Aizu Basin ranges from 1.5 × 10 −2 to 1.0 × 10 −1 m/d.The inclination of the slope in the Aizu Basin ranges from 10/1000 to 20/1000, except along the Aga River, which is located in central part of the basin [29].On the other hand, the average inclination of the slope in the Sendai Plain is 6/1000, and 0.4/1000 at alluvial fans existing in the inland area, and low-lying shore and dunes existing in the seacoast area [20].The smaller inclination of the slope of the plain may lead to lower groundwater flow velocity compared to that of the basin.The Péclet number (Pe) was calculated for each BHE location in the plain and basin to examine the difference of the ratio of advection by that of diffusion between the two.Pe is a dimensionless number and was calculated as follows: where ν is the average groundwater flow velocity at 0-100 m depth [LT −1 ], L is the representative length (=0.281 m, which is half the circumference of an actual GSHP system in Japan) [L], and α is average thermal diffusivity at 0-100 m depth.α was obtained by dividing the average thermal conductivity at 0-100 m depth by the volumetric heat capacity [L 2 T −1 ].For the BHEs located in the plain, Pe ranges from 3.5 × 10 −3 to 7.3 × 10 −2 , while the corresponding range for the basin is 8.5 × 10 −2 to 5.8 × 10 −1 .This indicates that compared to the Aizu Basin, heat transfer in the Sendai Plain is less affected by advection of groundwater flow.Therefore, the heat exchange rate may be impacted not only by groundwater flow velocity but also by subsurface temperature of the plain, and the heat exchange rate may be influenced only by groundwater flow velocity of the basin.It is unclear whether the heat exchange rate is affected by groundwater flow velocity and/or subsurface temperature (heat conduction) depending on Pe.This topic will be explored in a future study.The characteristics of the distribution of heat exchange rate were examined for the Sendai Plain and Aizu Basin.The upstream area is expected to have a higher heat exchange rate and would thus be more suitable for a GSHP system utilizing groundwater flow compared to the downstream area for both the plain and the basin.This is because the upstream area may have higher groundwater flow velocity than the downstream area.All basins (inland basins) do not necessarily have a higher inclination of slope than the plains (coastal plains), and actual groundwater flow velocity depends not only on the hydraulic gradient but also on hydraulic conductivity, which is controlled by geological facies [26].However, basins are expected to have higher groundwater flow velocity and heat exchange rate than plains in terms of the inclination of slope.This indicates that basins are expected to be more suitable for GSHP system installation due to the groundwater flow velocity.Seacoast areas of plains may have lower groundwater flow velocity because of the smaller inclination of slope, and utilizing aquifer thermal energy storage (ATES) systems is more suitable under such scenarios.

Conclusions
In this study, 3D groundwater flow and heat transport simulations were carried out for the Sendai Plain, the aim being to assess the installation potential of a GSHP system in the plain and compare the same with corresponding results of a past study on the Aizu Basin.The setting and results of this study are enumerated below: (1) The average heat exchange rate for the 33 constructed BHEs for the Sendai Plain was 22.5 W/m.
In this analysis, the operating scenario of the GSHP system was set as 120 d of space heating per year from December to March, assuming 24-h operation.The inlet temperature and flow rate of circulation fluid were set as 5 • C and 20 L/min, respectively.
(2) The Sendai Plain has lower heat exchange rates than the Aizu Basin.In both areas, upstream locations showed a higher heat exchange rate than downstream locations.(3) Multiple regression analysis was conducted using heat exchange rate as the response variable and average groundwater flow velocity and average subsurface temperature profile at 0-100 m depth as the explanatory variables.The heat exchange rate of the Sendai Plain, whose Péclet number ranges from 3.5 × 10 −3 to 7.3 × 10 −2 , may be affected by not only groundwater flow velocity but also by subsurface temperature.The heat exchange rate of the Aizu Basin, whose Péclet number ranges from 8.5 × 10 −2 to 5.8 × 10 −1 , may be affected by groundwater flow velocity alone.(4) In terms of inclination of slope, the upstream area is expected to have a higher heat exchange rate, which is more suitable for GSHP system installation, which can utilize the groundwater flow, compared to the downstream area in both the plain and the basin.Similarly, inland basins are expected to be more suitable for GSHP system installation utilizing groundwater flow velocity than the coastal plains.The area along the coast is expected to be more amenable to the use of ATES systems.
This study showed that multiple regression analysis can reveal factors affecting the heat exchange rate as well as the degree to which they affect it.Note, however, that this study is a comparative study, and a common conclusion cannot be drawn for other plains and basins due to their different geographical and hydrogeological conditions.However, this approach is novel in that it can be used to examine factors affecting the heat exchange rate in different regions of Japan and other countries.
This study conducted a suitability assessment of the GSHP system for space heating.It was found that subsurface temperature and outlet temperature of the GSHP system may influence the heat exchange rate in areas where this rate is affected by heat conduction.Further study in this field will focus on the installation potential of the system for space heating as well as space cooling in the same region, and the estimated distribution map of the heat exchange rate will be compared for these operations.

Figure 1 .
Figure 1.Study area and modeled area.The relief map and elevation data are drawn using the 10-m grid digital elevation map sourced from the Geospatial Information Authority of Japan (GSI).
the correction sink/source term of extended Oberbeck-Boussinesq approximation [T −1 ], K is the tensor of hydraulic conductivity [LT −1 ], f u is the viscosity relation function[-], χ e is the dimensionless buoyancy coefficient[-], is the porosity [-], ρ is the density of the fluid [ML −3 ], c is the specific heat capacity of the fluid [L 2 T −2 K −1 ], ρ s is the density of the solid [ML −3 ], c s is the specific heat capacity of the solid [L 2 T −2 K −1 ], T is the temperature [K], Λ is the tensor of hydrodynamic thermodispersion [MLT −3 K −1 ], Q T w is the specific heat sink/source function of wells [ML −1 T −3 ], T 0 is the reference temperature [K], Λ 0 is the isotropic thermal conductivity of the fluid [MLT −3 K −1 ], Λ s 0 is the isotropic thermal conductivity of the solid [MLT −3 K −1 ], D mech is the tensor of mechanical dispersion [L 2 T −1 ], Λ is the coefficient of thermalconductivity of the fluid [MLT −3 K −1 ], δ is the Dirac delta function [-], Λ s is the coefficient of thermal conductivity of the solid [MLT −3 K −1 ], β T is the transverse dispersivity [L], β L is the longitudinal dispersivity [L], H is the modified heat source/sink term without well function [L 2 T −3 ], and H s is the modified heat source/sink term of well function [L 2 T −3 ] [21].The terms listed on the right-hand side in Equations (1) and (3) are to consider groundwater discharge.

Figure 2 .
Figure 2. 3D groundwater flow and heat transport model of the Sendai Plain.

Table 1 .
Input parameters of the regional model of the Sendai Plain.

Figure 3 .
Figure 3. Location of groundwater monitoring wells in the Sendai Plain.

Figure 4 .
Figure 4. Comparison of measured hydraulic heads with calculated values.

Figure 5 .
Figure 5.Comparison of measured subsurface temperature profiles with calculated values.

Figure 6 .
Figure 6.Comparison of measured outlet temperature with calculated values.

Figure 7 .
Figure 7. Locations of the ground heat exchanger (GHE) models.

Figure 8 .
Figure 8. Distribution map of heat exchange rates for the Sendai Plain.

Figure 9 .
Figure 9. Distribution map of average groundwater flow velocity from the surface (0 m) to 100 m depth for the Sendai Plain.

Figure 10 .
Figure 10.Distribution map of average subsurface temperature from the surface (0 m) to 100 m depth for the Sendai Plain.

Figure 11 .
Figure 11.Distribution map of heat exchange rates in the Aizu Basin quoted from Shrestha et al. [19].

Figure 12 .
Figure 12. Results of the regression analysis.Relationship between average groundwater flow velocity and heat exchange rate for the (a) Sendai Plain and (b) Aizu Basin.Relationship between average subsurface temperature and heat exchange rate for the (c) Sendai Plain and (d) Aizu Basin.Results of the multiple regression analysis for the (e) Sendai Plain and (f) Aizu Basin.Data for the Aizu Basin were sourced from Shrestha et al. [16].

Table 2 .
Results of the multiple regression analysis.X 1 and X 2 denote average groundwater flow velocity and average subsurface temperature, respectively.