Proposal for a Method Predicting Suitable Areas for Installation of Ground-Source Heat Pump Systems Based on Response Surface Methodology

: The installation potential of ground-source heat pump (GSHP) systems has been studied based on the spatial interpolation of numerical simulation results using ground heat exchanger (GHE) models. This study is the ﬁrst to create an estimation formula for the heat exchange rate (HER) to obtain a solution equivalent to the numerical analysis results considering the average method when supplying three-dimensional (3D) hydrogeological information that a ﬀ ects the HER to a two-dimensional (2D) map. It was found that the main factors a ﬀ ecting the HER were groundwater ﬂow velocity, subsurface temperature, and thermal conductivity. The response surface methodology was utilized to approximate the HER using the above-mentioned three parameters. The estimated HER showed very strong agreement with that calculated by the GHE models. The application of the estimation formula to the simulation of the 3D groundwater ﬂow and heat transport model of the Sendai Plain (Japan) better reﬂects the hydrogeological information of the regional model than conventional maps. The proposed method improves the spatial resolution of maps and allows for the easy creation of the HER estimation formula.


Introduction
Ground source heat pump (GSHP) systems are an energy-efficient and environment-friendly technology [1][2][3][4]. They use the natural thermal energy of the subsurface stored at shallow depths for space heating/cooling and snow-melting, among other uses [5][6][7]. It is reported that the efficiency of thermal exchange of the closed loop system depends on thermal conductivity, heat capacity, and subsurface temperature profile [8]. Sensitivity analysis also showed the results that both the design parameters and operational parameters have a major influence on the GSHP performance [9]. The influence of unsaturated condition and groundwater flow on the performance of GSHPs was examined, which showed that groundwater level affected the performance of the GSHP system [10]. Therefore, it is important to subsurface conditions. In Europe, the "Cheap-GSHPs" (Cheap and efficient application of reliable Ground Source Heat exchangers and Pumps) European project has developed a Decision Support System to assist the decision-making process of GSHP system designers and building owners [11].
In Japan, GSHP suitability mapping has been carried out for increasing the number of installations of GSHP systems. The installation potential of GSHP systems was investigated using field survey data [12] and numerical simulation [1,13,14]. Many studies showed that groundwater flow velocity increases the heat exchange rate (HER) and temperature distribution in the vicinity of the ground heat exchanger (GHE) [10,[15][16][17][18][19]. The development and uptake of this system in Japan are gradually increasing. However, the number of installations of closed-loop GSHP systems in Japan is about 2700 as of March 2018 [20], which is lower than at other locations around the world, such as Europe and America. This is mainly because it is expensive to install such systems, and GSHP systems are not recognized in Japan. While policies on the introduction of renewable energies are being implemented in each municipality, GSHP systems have not been positioned as an environmental policy in many municipalities. It is also a problem that the public does not have the opportunity to know the GSHP systems. As one of the ways of solving this, suitability maps for the systems are created for each municipality or basin and plain in Japan. This map allows you to grasp the area of a region with relatively high potential for the GSHP systems, and it is expected that local governments and residents will be able to start and help when considering where to introduce the systems. In addition, it is also expected that publishing the map will lead the public to recognize the systems. This potential represents the result of HER under certain operating conditions of a GSHP system, and the method of evaluating this potential is as follows: first, a regional groundwater flow and heat transport simulation is performed. Next, the locations where the GHE model is created are selected in the area where the potential map is to be created, and the HER is calculated by inputting the simulation results and input parameters extracted from the regional model to each constructed GHE model. Finally, the potential map is created by spatially complementing the results of the GHE model. However, numerical simulation of the installation potential was done based on spatial interpolation of the HER in a geographic information system (GIS), and suitability maps for GSHP systems were created without fully reflecting the results of the numerical simulation. Increasing the number of pilot points increases the resolution of the potential map; however, it also increases the time needed to create the map, which is the challenge addressed in this study.
Here, a new method of evaluating the suitability of GSHP systems was proposed to reflect the hydrogeological conditions on the suitability map. Specifically, based on the numerical simulations which were conducted by Kaneko et al. [21], an optimization technique was developed and applied to the model area in this study. Kaneko et al. [21] constructed a regional groundwater flow and heat transport model, and heat exchange simulations have been conducted for each GHE model at the pilot points (GHE locations). We apply the mathematical optimization method to the numerical simulation results of the GHE locations to approximate the HER based on the numerical simulation results of the regional model. As a result, it is possible to create a map of the HER using the result of each mesh of the regional model with data on hydrogeological conditions. When creating a potential map using an optimization method, it is necessary to examine the average in the depth direction to transfer 3D hydrogeological data onto a 2D map, because this method uses average hydrogeological parameters. In this study, the hydrogeological parameters affecting the HER were first examined in terms of the average in vertical direction, which is the novelty of this study. Then, an estimation formula for the HER was created via the optimization method using the relationship between the HER calculated by the GHE models created by Kaneko et al. [21] and the hydrogeological parameters of the GHE models. These parameters were also extracted from the central coordinates of the regional model created by Kaneko et al. [21], and the estimation formula was applied to the extracted parameters to obtain a HER in the study area. The proposed suitable evaluation method, which combines numerical simulation and mathematical optimization methods, is unique to this study. Additionally, this method can be evaluated with higher accuracy and in a shorter time than before, which can be effective to select suitable locations for GSHP system installations and promote the growth of such systems in Japan.

Review of Hydrogeological Information and Numerical Simulation of the Study Area
The study area is the Sendai Plain, which is located in Miyagi Prefecture. The plain stretches in the N-S and W-E direction with a distance of 40 km and 10 km, respectively ( Figure 1). Small-scale  [21]. 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).  [21]. 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). Heat exchange simulations were conducted for space heating for each GHE model at the respective locations by Kaneko et al. [21] using FEFLOW [22]. Identical 20 m × 20 m × 120 m GHE models were constructed at 33 locations in the Sendai Plain, assuming a general closed-loop system ( Figure 3). Parameters input to the GHE models are the same as those assigned for the analysis model, as shown in Table 1. The same hydrogeological parameters from the same GHE locations in the regional model were assigned for each GHE model. Similarly, initial and boundary conditions of the GHE models were set based on the groundwater flow velocity, hydraulic head, and underground temperature calculated from the analysis model at each GHE location. A GHE was set at a depth of 100 m at the central nodes of the GHE models. A double U-tube with a diameter of 34 mm and silica sand grout were considered. The HER was calculated at each GHE location based on the following conditions: an operating scenario of 120 d of space heating per year from December to March, assuming 24-h operation; inlet temperature and flow rate of the circulation fluid at 5 °C and 20 L/min, respectively. In this paper, the operating scenario of the GSHP system is calculated by the above condition. This study was performed based on a regional model which has been calibrated, and the GHE models in the Sendai Plain (Japan) by Kaneko et al. [21]. The reginal model was verified by 1 TRT, 9 groundwater levels, and 11 subsurface temperature profiles. Heat exchange simulations were conducted for space heating for each GHE model at the respective locations by Kaneko et al. [21] using FEFLOW [22]. Identical 20 m × 20 m × 120 m GHE models were constructed at 33 locations in the Sendai Plain, assuming a general closed-loop system ( Figure 3). Parameters input to the GHE models are the same as those assigned for the analysis model, as shown in Table 1. The same hydrogeological parameters from the same GHE locations in the regional model were assigned for each GHE model. Similarly, initial and boundary conditions of the GHE models were set based on the groundwater flow velocity, hydraulic head, and underground temperature calculated from the analysis model at each GHE location. A GHE was set at a depth of 100 m at the central nodes of the GHE models. A double U-tube with a diameter of 34 mm and silica sand grout were considered. The HER was calculated at each GHE location based on the following conditions: an operating scenario of 120 d of space heating per year from December to March, assuming 24-h operation; inlet temperature and flow rate of the circulation fluid at 5 • C and 20 L/min, respectively. In this paper, the operating scenario of the GSHP system is calculated by the above condition. This study was performed based on a regional model which has been calibrated, and the GHE models in the Sendai Plain (Japan) by Kaneko et al. [21]. The reginal model was verified by 1 TRT, 9 groundwater levels, and 11 subsurface temperature profiles.   Figure 4 shows a flowchart of the proposed method in this paper. In this study, the hydrogeological condition input to the GHE model and the simulation result of the model were extracted, and the factors affecting HER were statistically analyzed. In the GHE model,   Figure 4 shows a flowchart of the proposed method in this paper. In this study, the hydrogeological condition input to the GHE model and the simulation result of the model were extracted, and the factors affecting HER were statistically analyzed. In the GHE model, hydrogeological parameters include horizontal/vertical hydraulic conductivity, specific storage, thermal conductivity of solid and fluid, volumetric heat capacity, and porosity. The simulation results also include hydraulic heat, subsurface temperature, and groundwater flow velocity. These parameters were average values from the surface to 100 m depth in the analysis because the length of the GHE was set to 100 m. As a result, the HER was correlated with groundwater flow velocity, v, subsurface temperature, T, and, thermal conductivity of the solid, λ s . It was judged that the HER can be expressed as the function which uses the above three variables. The response surface methodology (RSM) was utilized as a mathematical optimization method, which was used to approximate the HER by Tomigashi et al. [23], consisting of the above three variables in this study. When a quadratic polynomial is used as the response function, the response surface is expressed by the following equation:

Response Surface Methodology
where, β 0 , β i , β ii , and β ij are the regression coefficients and the x values are the variables. When the HER Q is approximated by the three variables, Equation (1) is expressed by the following biquadratic polynomial: where β 1 to β 9 are regression coefficients, v is the groundwater flow velocity, T is the subsurface temperature, and λ s is the thermal conductivity of solids, respectively. In this method, it is necessary to average the three parameters that change in the depth direction because the parameters input to the GHE model are not uniform values but values that can change in the depth direction. It is also reported by Vu-Bac et al. [24] that surrogate models such as polynomial regression models have to be utilized as alternative models based on which the sensitivity indices are estimated. In this study, the GHE model was created, and the averaging in the depth direction was examined before estimating the heat exchange amount using Equation (2) in the next section. This examination used the finite element software FEFLOW [22].

Groundwater Flow Velocity
The averaging of groundwater flow velocity was examined. The formations in the study area are mainly from the Quaternary and Neogene from the surface to a depth of 100 m, which were input to the regional model of the Sendai Plain by Kaneko et al. [21]. In this study, the averaging method assuming the above two layers was considered. The uniform model and two-layer model with the same arithmetic average groundwater flow velocity were constructed, and the HER was calculated. T was 15 • C and λ s was 1.5 W/(mK). Figure 5 shows the temperature distribution in the GHE of the uniform model, two-layer model, and HER when the arithmetic average groundwater flow velocity was 0.22 m/day constantly. The HER was 63.27 W/m when the v was 0.22 m/day (Figure 5a) constantly. When the velocity from the surface to 50 m depth (v 1 ) was 0.3 m/day and the velocity from a depth of 50 m to 100 m (v 2 ) was 0.14 m/day, the HER was 62.36 W/m ( Figure 5b). When v 1 was 0.4 m/day and v 2 was 0.04 m/day, the HER was 57.68 W/m ( Figure 5c). It was confirmed that even if the arithmetic average groundwater flow rate was the same, the temperature distribution of the inlet and outlet changed, and the heat exchange value was different. It was also found that the HER of the two-layer Energies 2020, 13, 1872 6 of 18 model was smaller than that of the uniform model, and the HER decreased as the difference between the v 1 and v 2 of the two-layer model increased. Figure 6 shows the relationship between groundwater flow velocity and HER in the uniform model. The relationship between the two was confirmed to be nonlinear. The HER of the two-layer model calculated under the conditions that v 1 was 0.3 m/day and v 2 was 0.14 m/day in Figure 5b corresponds to the HER of the uniform model with a groundwater flow velocity of 0.2 m/day. The HER of the two-layer model calculated under the conditions that v 1 was 0.4 m/day and v 2 was 0.04 m/day in Figure 5c also corresponds to the HER of the uniform model with the velocity of 0.13 m/day. This means that it is necessary to consider an averaging method that averages 0.3 m/day and 0.14 m/day to 0.20 m/day and that averages 0.4 m/day and 0.04 m/day to 0.13 m/day. To examine the above average method, the amount of heat exchange at that time was calculated by changing the v pattern of the uniform and the two-layer model. As a result, the average method can be explained as a combination of arithmetic average and harmonic average in the following equation: where v a is the arithmetic average groundwater flow velocity; v h is the harmonic average groundwater flow velocity; and α and β are parameters that vary with v a . The relationship between v a and the parameters α and β in Table 2 was interpolated with a cubic spline ( Figure 7). The average groundwater flow velocity of the two-layer model can be explained by an arithmetic average because the velocity is smaller; however, it is more affected by the harmonic mean when the velocity is greater.

Groundwater Flow Velocity
The averaging of groundwater flow velocity was examined. The formations in the study area are mainly from the Quaternary and Neogene from the surface to a depth of 100 m, which were input to the regional model of the Sendai Plain by Kaneko et al. [21]. In this study, the averaging method assuming the above two layers was considered. The uniform model and two-layer model with the same arithmetic average groundwater flow velocity were constructed, and the HER was calculated. T was 15 °C and λ s was 1.5 W/(mK). Figure 5 shows the temperature distribution in the GHE of the where is the arithmetic average groundwater flow velocity; is the harmonic average groundwater flow velocity; and and are parameters that vary with . The relationship between and the parameters and in Table 2 was interpolated with a cubic spline (Figure 7). The average groundwater flow velocity of the two-layer model can be explained by an arithmetic average because the velocity is smaller; however, it is more affected by the harmonic mean when the velocity is greater.

Subsurface Temperature
The average values of T were examined. Subsurface temperature is affected not only by heat dispersion but also by the advection of groundwater flow [25][26][27][28][29][30][31]. Subsurface temperatures in the geothermal zone, which occurs at depths between 10 and 20 m, normally follow the geothermal gradient which in turn is usually represented by an increase of 1 °C per 20 to 40 m [26,32,33]. Previous studies have noted temperature anomalies at certain depths due to the artificial pumping of groundwater and temperature inversions due to global warming [33]. This study considered the five subsurface temperature profile patterns shown in Figure 8 and investigated the effect that these patterns had on the HER. The arithmetic average value of all these profiles was 15 °C. v was 1e-6 m/day and λ s was 1.5 W/(mK). It was found that the HER did not change significantly in any profiles. Therefore, T can be calculated as an arithmetic average.

Subsurface Temperature
The average values of T were examined. Subsurface temperature is affected not only by heat dispersion but also by the advection of groundwater flow [25][26][27][28][29][30][31]. Subsurface temperatures in the geothermal zone, which occurs at depths between 10 and 20 m, normally follow the geothermal gradient which in turn is usually represented by an increase of 1 • C per 20 to 40 m [26,32,33]. Previous studies have noted temperature anomalies at certain depths due to the artificial pumping of groundwater and temperature inversions due to global warming [33]. This study considered the five subsurface temperature profile patterns shown in Figure 8 and investigated the effect that these patterns had on the HER. The arithmetic average value of all these profiles was 15 • C. v was 1e-6 m/day and λ s was 1.5 W/(mK). It was found that the HER did not change significantly in any profiles. Therefore, T can be calculated as an arithmetic average. studies have noted temperature anomalies at certain depths due to the artificial pumping of groundwater and temperature inversions due to global warming [33]. This study considered the five subsurface temperature profile patterns shown in Figure 8 and investigated the effect that these patterns had on the HER. The arithmetic average value of all these profiles was 15 °C. v was 1e-6 m/day and λ s was 1.5 W/(mK). It was found that the HER did not change significantly in any profiles. Therefore, T can be calculated as an arithmetic average.

Thermal Conductivity
The average of λ s was examined. The typical value of the unconsolidated sediment is from 1 to 2 W/(mK) [34]. Assuming the thermal conductivity shown in Figure 9, the effect of the change in the thermal conductivity in the depth direction on the HER was examined. The arithmetic average values of λ s were 1.5 W/(mK) in all scenarios, in which v was 1.0 × 10 −6 m/day and T was 15 • C. In this calculation pattern, a maximum difference of 1.6% occurred comparing Figure 9a,b. The thermal conductivity of the wide area model created by Kaneko et al. [21] was 1.4 W/(mK) for the Quaternary and 1.5 W/(mK) for the Neogene. Even if arithmetic averaging is performed, it is considered to have almost no effect on the HER. Therefore, λ s can be calculated as an arithmetic average.

Creation of Estimation Formula for the HER
An estimation formula for the HER was created using the HER calculated for 33 GHE models by Kaneko et al. [21], and the three parameters v, T, and λ s input to the GHE model in Sendai Plain. These parameters are average values calculated by the method shown in Section 3.2. In this study, the regression coefficients β 0 to β 9 were determined by the least squares method to create an approximate expression using Equation (2). Table 3 shows the parameters of the estimation formula, and Table 4 also shows relational data between the three parameters and HER at 33 GHE models. The HER was estimated by substituting the three parameters and regression coefficients into Equation (2). Figure 10 shows a comparison of the HER calculated by the GHE models, with one estimated by the biquadratic polynomial of the RSM (Equation (2)). The coefficient of determination between the HER calculated by the GHE models and those estimated by the RSM in the Sendai Plain was 0.999, and the root mean square error (RMSE) between these two was 0.129. It was judged that the HER estimated by the RSM was reasonable.
Energies 2020, 13, 1872 10 of 18 2 W/(mK) [34]. Assuming the thermal conductivity shown in Figure 9, the effect of the change in the thermal conductivity in the depth direction on the HER was examined. The arithmetic average values of λ s were 1.5 W/(mK) in all scenarios, in which v was 1.0 × 10 −6 m/day and T was 15 °C. In this calculation pattern, a maximum difference of 1.6% occurred comparing Figure 9a and 9b. The thermal conductivity of the wide area model created by Kaneko et al. [21] was 1.4 W/(mK) for the Quaternary and 1.5 W/(mK) for the Neogene. Even if arithmetic averaging is performed, it is considered to have almost no effect on the HER. Therefore, λ s can be calculated as an arithmetic average.

Creation of Estimation Formula for the HER
An estimation formula for the HER was created using the HER calculated for 33 GHE models by Kaneko et al. [21], and the three parameters v, T, and λ s input to the GHE model in Sendai Plain. These parameters are average values calculated by the method shown in Section 3.2. In this study, the regression coefficients to were determined by the least squares method to create an approximate expression using Equation (2). Table 3 shows the parameters of the estimation formula, and Table 4 also shows relational data between the three parameters and HER at 33 GHE models. The HER was estimated by substituting the three parameters and regression coefficients into Equation (2). Figure 10 shows a comparison of the HER calculated by the GHE models, with one estimated by the biquadratic polynomial of the RSM (Equation (2)). The coefficient of determination  between the HER calculated by the GHE models and those estimated by the RSM in the Sendai Plain was 0.999, and the root mean square error (RMSE) between these two was 0.129. It was judged that the HER estimated by the RSM was reasonable.    Figure 10. Comparison of the heat exchange rate calculated by the GHE models with that estimated by the RSM in the Sendai Plain. HER, GHE, and RSM denote the heat exchange rate, ground heat exchanger, and response surface methodology, respectively.

Application of the Estimation Formula for the HER to the Sendai Plain
The HER in the Sendai Plain was estimated based on the regional model by Kaneko et al. [21]. The average value of three parameters (v, T, and λ s ) from the surface to 100 m depth was extracted from meshes of the regional model, and the HER was estimated by substituting the three parameters into Equation (2). A distribution map of HER was created by performing spatial interpolation using the Kriging method on the HER estimated by each model mesh. Figure 11 shows a comparison of the distribution map of HERs by Kaneko et al. [21] (Figure 11a, conventional map) with that estimated by the RSM (Figure 11b, refined map) in the Sendai Plain. Figure 12 also shows distribution maps of average groundwater flow velocity (Figure 12a), subsurface temperature (Figure 12b), and thermal conductivity of solids (Figure 12c) from the surface (0 m) to 100 m depth for the plain. There is little difference in the HER between the two maps at the GHE locations; however, the tendency of distribution of the HER between GHE locations was different. This is because the conventional map was created by spatial interpolation based on the HER calculated by 33 GHE locations, whereas the refined map reflects the distribution of the three parameters of v, T, and λ s . It was found that the area where v or T was greater had a greater HER in the refined map. It was reported that the HER was strongly related to v and T for the plain by Kaneko et al. [21].

Validation of the Estimation Formula by Constructing GHE Models
To verify the validity of the refined map, GHE models were created at eight points (eight meshes of the regional model) where the difference in the HER between the conventional map and the improved map was relatively larger (Figure 13). The HER was calculated under the same operating conditions of the GSHP system as that used by Kaneko et al. [21]. Figure 14 shows a comparison of the HER calculated by the GHE model and the HER estimated by the RSM. The RMSE was 0.90 for the refined map and 4.30 for the conventional map. The coefficient of determination also was 0.988 for the refined map and 0.401 for the conventional map. This is because the conventional map is created by spatial interpolation based on the HER calculated at 33 GHE locations, while the refined map is created by spatial interpolation based on the HER estimated at 1273 meshes of the reginal model. Therefore, the refined map can better reflect the hydrogeological information of the regional model than the conventional map. This method cannot only create the estimation formula for the HER equivalent to the numerical analysis results but also increase the spatial resolution of the HER map.

Validation of the Estimation Formula by Constructing GHE Models
To verify the validity of the refined map, GHE models were created at eight points (eight meshes of the regional model) where the difference in the HER between the conventional map and the improved map was relatively larger (Figure 13). The HER was calculated under the same operating conditions of the GSHP system as that used by Kaneko et al. [21]. Figure 14 shows a comparison of the HER calculated by the GHE model and the HER estimated by the RSM. The RMSE was 0.90 for the refined map and 4.30 for the conventional map. The coefficient of determination also was 0.988 for the refined map and 0.401 for the conventional map. This is because the conventional map is created by spatial interpolation based on the HER calculated at 33 GHE locations, while the refined map is created by spatial interpolation based on the HER estimated at 1273 meshes of the reginal model. Therefore, the refined map can better reflect the hydrogeological information of the regional model than the conventional map. This method cannot only create the estimation formula for the HER equivalent to the numerical analysis results but also increase the spatial resolution of the HER map.

Discussion
By response surface methodology, which is one kind of optimization method and a surrogate model, the HER can be estimated without heat exchange simulation by GHE model calculation ( Figure 10). This proposed method can highly increase the spatial resolution of the HER map as compared to the conventional map ( Figure 11). In previous studies, increasing the number of pilot Figure 14. Comparison of the heat exchange rate calculated by GHE models with that estimated by the RSM in the Sendai Plain.

Discussion
By response surface methodology, which is one kind of optimization method and a surrogate model, the HER can be estimated without heat exchange simulation by GHE model calculation ( Figure 10). This proposed method can highly increase the spatial resolution of the HER map as compared to the conventional map ( Figure 11). In previous studies, increasing the number of pilot points increased the time needed to create the map. However, the proposed method cannot only estimate the HER from hydrogeological parameters in an easy way but also increase the spatial resolution of the map without performing heat exchange simulations by creating a GHE model. By using this method, it is proposed that the HER can be estimated in other regions if the three parameters of v, T, and λ s can be obtained. By using the proposed method, users will have the benefit of estimating the HER with field survey data without performing numerical simulations in the future.
However, there are some challenges for applying the method to other regions. In this study, regression coefficients were calculated using the relationship between the HER and the three parameters of the GHE model in the Sendai Plain. Thus, if the three parameters outside the range of the relational data were input, misestimation of the HER may occur. Focusing on the three parameters of the 33 GHE models used to create the regression coefficient of the estimation formula for the HER, v ranged from 9.5 × 10 −4 to 2.0 × 10 −2 m/day, T ranged from 13.1 to 20.8 • C, and λ s ranged from 1.4 to 1.5 W/(mK). When using the parameters of the regression coefficients in Table 3, the estimation formula can only be applied within the above ranges of the three hydrogeological parameters. For example, if v is set to 1.0 m/day, T is set to 15 • C, and λ s is set to 1.5 W/m/K, and the parameters and the regression coefficients of Table 1 are input to Equation (2), the HER is calculated as −14,085 W/m, which is an obvious mistake. As a countermeasure to this, increasing the relational data between the HER and the three variables is expected to create an estimation formula for the HER applicable to any region.
Moreover, what needs to be considered is the averaging method of the three-dimensional hydrogeological information in the depth direction, especially the averaging method of the groundwater flow velocity. In the Sendai Plain, formations in the area are mainly from the Quaternary and Neogene from the surface to a depth of 100 m, which were input to the reginal model of the Sendai Plain by Kaneko et al. [21]. In this study, the averaging method assuming the above two layers was considered. Under this condition, the HER can be estimated by using the averaging method of v, which is a combination of the arithmetic average and harmonic average shown in Equations (3) and (4).
Here, the proposed method was applied to 20 GHE models in the Aizu Basin, Japan constructed by Shrestha et al. [35]. The regression coefficients β 0 to β 9 were determined by the least squares method to create an approximate expression using Equation (2). Table 5 shows the parameters of the estimation formula. The HER was estimated by substituting three parameters and regression coefficients into Equation (2). Figure 15 shows a comparison of the HER calculated by the GHE models with that estimated by the biquadratic polynomial of the RSM (Equation (2)). The coefficient of determination between the HER calculated by the GHE models and those estimated by the RSM in the Sendai Plain was 0.936, and the RMSE between the two in the Sendai Plain was 1.100. The estimation accuracy of the HER in the basin was not better than that in the plain. This is because the input parameters to the GHE models in the vertical direction between the two regions were different. In the plain, the input parameters were changed by the boundary between Quaternary and Neogene. In contrast, they were changed by 10 m depth in the basin.   (Figure 16a) constantly. When the velocity from the surface to 50 m depth (v 1 ) was 0.4 m/day and the velocity from a depth of 50 m to 100 m (v 2 ) was 0.04 m/day, the HER was 57.68 W/m (Figure 16b). However, when v 1 was 0.4 m/day and v 2 was 0.04 m/day, which were entered alternately every 10 m in depth, the HER was 61.34 W/m (Figure 16c). The averages in the conditions of Figure 16b,c are both calculated as arithmetic average of 0.22 m/day and the harmonic average of 7.3 × 10 −2 m/day. The same HER was estimated using the optimization method under the conditions in Figure 16b,c. This is because Equations (3) and (4) for averaging v are calculated by a combination of arithmetic average and harmonic average; therefore, they have the same average groundwater flow velocity, that is, the same HER is obtained under the conditions in Figure 16b,c. Equations (3) and (4) are intended for the two-layer model (Sendai Plain hydrogeological model), so it is necessary to create another averaging equation not only for the 10-layer model but also for other models, which is a content of further study. Here, the proposed method was applied to 20 GHE models in the Aizu Basin, Japan constructed by Shrestha et al. [35]. The regression coefficients to were determined by the least squares method to create an approximate expression using Equation (2). Table 5 shows the parameters of the estimation formula. The HER was estimated by substituting three parameters and regression coefficients into Equation (2). Figure 15 shows a comparison of the HER calculated by the GHE models with that estimated by the biquadratic polynomial of the RSM (Equation (2)). The coefficient of determination between the HER calculated by the GHE models and those estimated by the RSM in the Sendai Plain was 0.936, and the RMSE between the two in the Sendai Plain was 1.100. The estimation accuracy of the HER in the basin was not better than that in the plain. This is because the input parameters to the GHE models in the vertical direction between the two regions were different. In the plain, the input parameters were changed by the boundary between Quaternary and Neogene. In contrast, they were changed by 10 m depth in the basin. Figure 15. Comparison of the heat exchange rate calculated by GHE models with that estimated by the RSM in the Aizu Basin. HER, GHE, and RSM denote heat exchange rate, ground heat exchanger and response surface methodology, respectively.   (3) and (4) for averaging v are calculated by a combination of arithmetic average and harmonic average; therefore, they have the same average groundwater flow velocity, that is, the same HER is obtained under the conditions in Figure 16b,c. Equations (3) and (4) are intended for the two-layer model (Sendai Plain hydrogeological model), so it is necessary to create another averaging equation not only for the 10-layer model but also for other models, which is a content of further study. Figure 15. Comparison of the heat exchange rate calculated by GHE models with that estimated by the RSM in the Aizu Basin. HER, GHE, and RSM denote heat exchange rate, ground heat exchanger and response surface methodology, respectively. From the above, using the proposed optimization method cannot only estimate the HER from hydrogeological parameters in an easy way but also highly increase the spatial resolution of a map by creating a GHE model without performing heat exchange simulation. The creation of the potential maps in various regions can be expected by using this method and by increasing the relational data between the HER and hydrogeological parameters. To create a universal estimation formula for the HER, it is necessary to consider the averaging method that matches the hydrogeological model of each region, which is a content of further study. From the above, using the proposed optimization method cannot only estimate the HER from hydrogeological parameters in an easy way but also highly increase the spatial resolution of a map by creating a GHE model without performing heat exchange simulation. The creation of the potential maps in various regions can be expected by using this method and by increasing the relational data between the HER and hydrogeological parameters. To create a universal estimation formula for the HER, it is necessary to consider the averaging method that matches the hydrogeological model of each region, which is a content of further study.

Conclusions
This study proposed an estimation formula for the HER to obtain a solution equivalent to the numerical analysis result considering the average method when supplying 3D hydrogeological information affecting HER to a 2D map. The results of the study are enumerated below: (1) It was found that the main factors affecting the HER were groundwater flow velocity (v), subsurface temperature (T), and thermal conductivity of solids (λ s ). (2) The average method in the vertical direction (GHE length) was evaluated. T and λ s can be calculated by arithmetic averaging. In contrast, v has to be calculated as a combination of arithmetic and harmonic average. (3) The RSM was utilized to approximate the HER using the three parameters in the Sendai Plain. The estimated HER agreed well with the calculated one by the GHE models because the coefficient of determination and RMSE between the two HERs were 0.999 and 0.129, respectively. (4) The proposed method cannot only estimate the HER from hydrogeological parameters in an easy way but also highly increase the spatial resolution of a map by creating a GHE model without performing heat exchange simulations.
It should be noted that the proposed method of suitable evaluation combining numerical simulation and mathematical optimization methods is unique to this study. This method can achieve higher accuracies in a shorter time compared to previous methods, which can help selecting suitable locations for GSHP system installations more effectively and promoting such systems in Japan.
For future study, an averaging method that matches the hydrogeological model of each region will be examined to create a universal estimation formula for the HER, leading to a potential map of GSHP systems throughout Japan that will increase the number of GSHP system installations.
Author Contributions: Conceptualization, S.K. and A.T.; methodology, S.K. and A.T. and M.Y.; investigation, T.I. and G.S.; writing-original draft preparation, S.K.; writing-review and editing, A.T.; visualization, S.K.; supervision, Y.U. and A.T. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.