Evaluating Variability of Ground Thermal Conductivity within a Steep Site by History Matching Underground Distributed Temperatures from Thermal Response Tests

Abstract

The variability of ground thermal conductivity, based on underground conditions, is often ignored during the design of ground-source heat pump systems. This study shows a field evidence of such site-scale variations through thermal response tests in eight borehole heat exchangers aligned at a site on a terrace along the foothills of mountains in northern Japan. Conventional analysis of the overall ground thermal conductivity along the total installation length finds that the value at one borehole heat exchanger is 2.5 times that at the other seven boreholes. History matching analysis of underground distributed temperature measurements generates vertical partial ground thermal conductivity data for four depth layers. Based on the moving line heat source theory, the partial values are generally within a narrow range expected for gravel deposits. Darcy velocities of groundwater are estimated to be 74–204 m/y at the borehole with high conductivity, increasing in the shallow layers above a depth of 41 m. In contrast, the velocities at the other seven boreholes are one-to-two orders of magnitude smaller with no trend. These high and low velocity values are considered for the topography and permeability. However, the relatively slow groundwater velocities might not apparently increase the partial conductivity.

1. Introduction

The ground thermal conductivity is a key determinant of the required size of a borehole heat exchanger (BHE) when designing and planning ground-source heat pump (GSHP) systems. This geothermal property is defined as a thickness-weighted average of the effective thermal conductivity of the soils/rocks comprising the formation at a site. The effective thermal conductivity of individual components ranges between <1 W/(m⋅K) for unsaturated fine sediments and >4 W/(M·K) for crystallized rocks, depending on the particles, porosity, and water content [1,2]. The design guideline of ASHRAE (American Society of Heating, Refrigerating and Air-Conditioning Engineers) specifies that a required length of a BHE in limestone with thermal conductivity of 3.1 W/(m·K) should be approximately half of that in unconsolidated clay of 1.2 W/(m·K) for the same conditions of thermal pulses in six hours [3]. The differences in BHE lengths became longer as the spans for thermal pulses became larger. Although GSHPs are established renewable-energy-based systems contributing to energy-saving and CO₂ emission reduction [4], the problem is that the capital costs of installing BHEs often exceed those of competing systems conventional air-source heat pumps [5]. Therefore, it is critical to obtain realistic estimates of the ground thermal conductivity at a site to determine the required lengths of BHEs, especially in the initial stage of designing and planning a heating, ventilation, and air-conditioning system [6].

Thermal response testing (TRT) is a standard method for evaluating ground thermal conductivity in-situ. It was originally developed in Sweden [7], and has since been standardized in various countries, such as the USA [8], Germany [9], and Japan [10], and by the international association [11]. The measurement and analysis methodologies are simple for practice. The test is conducted by circulating fluids in BHEs at a constant heating rate, and the temperature of the circulating fluids is observed at the shanks of inlet and outlet U-tubes. The temperature changes were analyzed based on heat-transfer theory assuming that the testing BHE is an infinite line source (ILS) in the ground [12]. This conventional TRT can provide ground thermal conductivity straightforwardly without significant computational effort.

In several decades, the methodology has been developed to improve its accuracy by incorporating other theoretical solutions for infinite cylindrical sources [1,11,13], and finite line sources [11,14]. The methodology used for determining ground thermal conductivity is an uncertain variable that could be modified to obtain sufficient agreements between measurements and calculations. The sum of squared error or the root mean squared error in the calculations is usually used for evaluating the agreements [15]. If a temporal series of fluid and/or soil temperatures are compared between calculations and measurements, the inversion method is called history matching analysis [13]. The literature review found that the share of applications of cylindrical and line source theory were 25% and 75%, respectively [16]. Complicated test conditions have also been accommodated by adopting numerical simulation [17,18,19]. The numerical modeling is available to not only estimate ground thermal conductivity, but also other parameters, such as thermal capacity [20], and borehole thermal resistance [21]. The numerical models are also used to estimate spatial distributions of ground thermal conductivity in a region [22]. Another testing approach based on variable heating power was also proposed in the same concept of water well pumping tests [23].

Ground thermal conductivity is potentially variable depending on the geological and groundwater conditions. A known issue in shallow geothermal utilization is uncertainty about such underground variation within site. TRT provides only one value at a location for the installation length of a test BHE, and only one test is usually conducted for cost reasons. Given the limited information, site-scale variability is assumed to be negligible in the practical design and planning of GSHP systems. This assumption might be valid under the following three conditions: The installation length of the BHE is fixed to neglect vertical changes; the formation is stratified, and each layer of a particular geologic component is continuous for negligible changes in horizonal directions; moreover, groundwater flows are slow enough to neglect advective heat transfer in the soil and rock.

The appropriate length of a BHE can be found for a site-dependent condition of vertical change in ground thermal conductivity through analysis of underground data collected by distributed temperature sensors (DTSs). DTSs are installed in a BHE to measure continuous depth-resolved temperature measurements during the heating phase, as called distributed TRT (DTRT). DTRT is classified in terms of installation number and location of DTSs; (1) a single DTS outside U-tubes [24], (2) a looped DTS attached along U-tubes [25], (3) a DTS into an inlet or outlet U-tube [13], and (4) a pair of DTS into inlet and outlet U-tubes [26,27]. Instead of circulating fluids at constant heating, the electric current-carrying copper wire is also used as a heat source along the DTSs [28,29,30]. The studies above have demonstrated to estimate vertical profiles of ground thermal conductivity from the distributed temperature data. A part of the studies also analyzed the profiles of borehole thermal resistance [26]. The heating cable TRTs were available in the boreholes at geotechnical investigations to cut the drilling costs [31].

This study focuses on the advection in TRTs by groundwater flows, the third listed effect on ground thermal conductivity. Previous field and laboratory studies of TRTs have shown advection to increase the ground thermal conductivity [32]. As an example, the apparent value was 9.14 W/(m·K) at a site of high groundwater flows in Japan [33]. The site was located in a basin surrounding the mountains with large precipitations, and natural groundwater flows at a site were estimated at 2750 m/y in the shallow layer [34]. Artificial flows of groundwater occur around water wells, potentially increasing the ground thermal conductivity. The authors demonstrated the DTRT near the pumping well, indicating the value increased from 2.00 W/(m.K) to 3.19 W/(m·K) [35]. In the guideline [36], the specific heat extraction per meter in depth is 80–100 W/m for “strong groundwater flow in gravel and sand,” and the value is almost double the value for no flow, 55–65 W/m. This means that the required length of a BHE is potentially designed to be 200% shorter, if groundwater flows actively. A theoretical study [37] indicated that the advection effect is not negligible when the Darcy velocity (specific discharge) is over 10⁻⁸ m/s (2 m/y), which is often the case in sedimentary basins. However, an actual effect of groundwater flows on the life cycle costs of a household GSHP system in Japan would appear when the Darcy velocity of groundwater is the order of 10¹ m/y or more [38]. The advection effects were also validated for multiple BHEs in numerical studies [39,40].

Although several research papers regarding TRTs was available, few studies have examined the variability of ground thermal conductivity on the intra-site scale for the construction of multiple BHEs, i.e., with an interval of several meters in a range of 10¹ to 10² m. One reason for the lack of relevant field data is the additional cost of conducting multiple TRTs at each site, as described above. Furthermore, the site topography might limit studies, with building sites often being located on flat land with slopes too gentle for significant advection. For example, the previous study demonstrated TRTs at six BHEs of different configurations in terms of U-tube and borehole diameter, indicating the small variability of ground thermal conductivity between 1.2–1.6 W/(m·K) [41]. However, in the steep reliefs, as often seen in Japan and other mountainous and hilly countries, high hydraulic gradients make rapid groundwater flows common, as described above. In such a site of high groundwater flows, it is expected to design reasonable lengths of BHEs, compared with no groundwater flows [33,34]. On the other hand, we also assumed that the groundwater flows might be uniform in space depending on the porosity distribution in the aquifer, because the sediments become coarser and more permeable around the foothills of mountains [42]. In such uniformity of groundwater flows, ground thermal conductivity is apparently changed. Consequently, the advection might affect more significantly the performances of multiple BHEs than in the flat plains.

This study aimed to evaluate site-scale variability of the ground thermal conductivity through TRTs at eight BHEs aligned at a site in Bibai City, Japan. The site was located on the terraced foothills of nearby mountains. The original purpose of the DTRTs was to compare borehole thermal resistances with respect to thermal performances of different configurations of U-tube, as previously reported [43]. The previous study summarized the ground thermal conductivities of the eight BHEs, indicating that the value at one site was about 2.5 times larger than those at the other seven BHEs. This study demonstrated a method to estimate the vertical changes of partial ground thermal conductivity in each of four layers comprising the formation through history matching analysis of underground distributed temperature data from optic-fiber DTSs. The present study compared the theoretical solutions for an ILS and a moving line source (MLS) [44]. The latter solution provided the partial ground thermal conductivity and Darcy velocity in each layer of the eight BHEs. The statistics represent field evidence for the site-scale variability in properties around steep reliefs.

2. Site Descriptions and Test Conditions

The study site was located at 43.32° N, 141.85° E, in Bibai City on Hokkaido, the northernmost island of Japan, as shown in Figure 1a. Figure 1b shows the Ishikari River, one of the island’s longest rivers, flowing from NNE to SSW from the northern mountains. The river’s valley has braided zones and alluvial fans consisting of coarse sediments in the Quaternary and Tertiary periods [45]. Near Sapporo, the island’s capital (population, about 2 million), the river’s flow changes to NE, before reaching the Japan Sea. There are meandering zones of Holocene fine sediments in the fluvial areas around Sapporo.

The study site is located on a small terrace along the foothills of the mountains. Drilling report of the BHEs construction showed that the formation was composed of monotonic gravel deposits, which were transferred mainly from the upper mountains, with some coming from the Ishikari River. Although in-situ measurements of permeability around the site were not obtained, the coarse deposits were assumed to be highly permeable. Also, the sediments were assumed not to be sorted well in the steep topography. This means the porosity related to soil particle distributions might be strongly uniform in the gravel deposits.

Eight BHEs, Nos 1 to 8, were constructed at the site in September 2018 at regular intervals of 4.0 m (except for the 6.5 m distance between Nos 6 and 7), as shown in Figure 2a. Each commonly had an installation depth of 82 m. The U-tubes were made of high-density polyethylene (HDPE), and were buried with silica sand, a typical grout material in Japan. The sand’s thermal properties were assumed to be close to those of the underground gravel deposits. The borehole diameter was 0.126 m for the single U-tubes of Nos. 3 to 8, and 0.139 m for double U-tubes. The eight BHEs were originally constructed to compare different cross-sectional configurations of U-tubes: Conventional circular shape or oval shape, single or double, and with or without a spacer to prevent thermal interference between the inlet and outlet U-tubes. The spacer length Ls was changed to evaluate its effect on thermal performance. Figure 2b shows the schematic cross-section of each BHE: No. 1, circle/double (Cd); No. 2, Cd with a spacer (82 m long) (Cd+); No. 3, circle/single (C); No. 4, C with a spacer (C+); No. 5, oval/single (O); No. 6, O with a spacer (Ls = 10 m, O+); No. 7, O+ (40 m); No. 8, O+ (82 m). The oval U-tube (O) was developed to improve efficiency by increasing the tube outer surface areas in boreholes with the same diameter as Cs boreholes [46].

Optic-fiber DTSs were attached to the surface of inlet and outlet U-tubes installed in each borehole. Each sensor was made of GI62.5 optic fiber (diameter, 125 μm) covered by a spiral metal tube for protection against soil pressure and by flexible plastic for water resistance, as shown in the upper right of Figure 2b. As the sensors were each attached twice, the sensor length in each BHE was 164 m (82 × 2). Each of the two attachments was cross-sectionally opposite. After construction of all the BHEs, the fibers were connected as a sequence from No. 1 to 8, resulting in a total length of over 1300 m. The temperatures on the outer surface of each U-tube were measured at regular 0.5 m depth intervals by the DTSs. The temperatures on the surface were assumed to be those at two opposite points at the same distance from the center, r_o. The DTS measurements were performed by using the signal receivers (N4385A type of A.P. Sensing, Co., Ltd., Boeblingen, Germany). The equipment has the advantage of changing the direction to send the optical signals in the DTS regularly to eliminate the fluctuation errors. The errors were also reduced through the filtering based on the Fast Fourier Transform (FFT) procedure, resulting in an averaged accuracy of 0.25 K, as shown in the previous study [47]. This study’s history matching analysis used the average of two temperatures by the DTSs, recorded on the opposing points at each depth T₀.

This study conducted conventional TRTs at each of the eight BHEs, as shown in Figure 3. Tap water was injected to fill the U-tubes at a pressure of about 0.1 MPa. The TRT equipment consisted of an electric heater, a circulation pump, a magnetic flow meter, and data acquisition devices. The water was circulated constantly at a flow rate of about 20 L/min. This study used two sets of equipment to save time. The circulating water was also heated at a constant rate of about 4.0 kW in one set of TRT equipment or 4.4 kW in the other. The heating period was 48 h (2 days), except at BHE No. 3 for high ground thermal conductivity, whose period was set at 96 h (4 days) to obtain stable temperature changes, likely due to groundwater flows. The temperatures of the circulating fluids were observed every minute at the shanks of the inlet and outlet U-tubes by Pt100 temperature sensors with an accuracy of 0.1 K.

The testing time schedule was as follows: 2018-10-10 16:01 to 2018-10-12 16:01 in Nos 1 and 3, 2018-10-23 12:31 to 2018-10-25 12:31 in Nos 2 and 6, 2018-10-16 15:31 to 2018-10-18 15:34 in Nos 4 and 8, 2018-10-31 12:53 to 2018-11-2 12:34 in No. 5, and 2018-11-6 15:02 to 2018-11-8 15:19 in No. 7. After heating, observation by the DTSs was maintained to measure recovery to the initial temperature.

3. Methodology for Analysis

First, this study performed conventional TRT analysis to determine overall values of ground thermal conductivity for the total lengths (82 m), $\overline{\lambda }$ in the eight BHEs. The analysis assumed each BHE to be an ILS with a constant heat exchange rate q [W/m], which was calculated from the total input power $Q$ [kW] divided by the installation length L [m]. The total input power $Q$ was calculated from the fluid temperature difference between the inlet and outlet shanks:

$$Q={\rho }_{\mathrm{w}}{C}_{\mathrm{w}}v\left(T_{f-in}-T_{f-out}\right),$$

(1)

where ${\rho }_{\mathrm{w}}{C}_{\mathrm{w}}$ is the volumetric heat capacity of the circulating water [J/(m³⋅K)]. $v$ is the circulating flow discharge [m³/s]. The input power $Q$ was confirmed to be kept within a 5% range, due to the fluctuations of $T_{f-in}$$T_{f-out}$, and $v$, and thus, the average Q was used to calculate q for analysis. We also assumed that the surrounding soil was homogeneous and isotropic. Thus, the change in underground soil temperature $\Delta T_s$ [K] by heating was calculated as [1]

$$\Delta T_s=T_s-T_0=\frac{q}{4\mathsf{\pi }\lambda }{{\displaystyle \int }}_{\frac{r^2}{4at}}^{\infty }\raisebox{1ex}{$e^{-u}$}\!\left/ \!\raisebox{-1ex}{$u$}\right.du,$$

(2)

where $T_s$ is the borehole surface temperature [°C], $T_0$ is the initial soil temperature [°C], $\lambda$ is the ground thermal conductivity [W/(m⋅K)], $a$ is the thermal diffusivity [m²/s] equal to $\lambda$ divided by volumetric thermal capacity—which is assumed to be 3 MJ/(K·m³)—$r$ is radius [m], and $t$ is elapsed testing time [h]. After the time greater than the minimum time, e.g., $t$ > 20 $a$/$r^2$ [1], Equation (2) could be approximated for simple calculation as

$$\Delta T_s\approx \frac{q}{4\mathsf{\pi }\lambda }\left[\log \left(\frac{4at}{r^2}\right)-\mathsf{\gamma }\right],$$

(3)

where $\mathsf{\gamma }$ is the Euler constant. Equation (3) indicates that the temperature increased log-linearly with time. Thus, the linear slope $\kappa$ was determined from least-squares fitting of semi-log plots of the average inlet and outlet temperatures from the PT100 sensors, i.e., $\overline{T_f}=\left(T_{f-in}+T_{f-out}\right)/2$. Considering the coefficient $q/4\mathsf{\pi }\lambda$ in Equation (2) was equal to the slope $\kappa$ of the measured logarithmic temperatures, $\overline{\lambda }$ was determined as

$$\overline{\lambda }=\frac{q}{4\mathsf{\pi }\kappa },$$

(4)

each BHE was not actually an ILS, and the temperature change of the circulating fluid $\overline{T_f}$ was affected not only by $\overline{\lambda }$, but also by the borehole thermal resistance R_b [m·K/W]. The value of $R_{\mathrm{b}}$ was calculated for the temperature difference between the fluid and borehole surface as

$$R_{\mathrm{b}}=\frac{T\prime -T_0}{q}-\frac{1}{4\mathsf{\pi }\overline{\lambda }}\left[\ln \frac{4a\times 3600}{r_{\mathrm{b}}{}^2}-\mathsf{\gamma }\right],$$

(5)

where $T\prime$ is the intercept of the fitted straight line [°C], for determination of $\kappa$ above, and represents the temperature at $t$ = 1 h ($\log t$ = 0). The blank value for 3600 was used for an elapsed time of 1 h (3600 s).

Estimation by history matching analysis was next applied using the distributed temperature measurements on the outer surfaces of the U-tubes during heating and recovery. The analysis obtained the vertically partial ground thermal conductivity ${\lambda }_{\mathrm{I}}, {\lambda }_{\mathrm{II}}, {\lambda }_{\mathrm{III}}, \mathrm{and} {\lambda }_{\mathrm{IV}}$ for layers I (depth; Z = 0–20.5 m), II (20.5–41 m), III (41–61.5 m), and IV (61.5–82 m), respectively. Each layer was assumed to be of equal thickness owing to limited information about the geologic structure, as the monotonic gravelly formation. The heat exchange rate q was also estimated to be constant within, but different among the layers. Vertically partial estimates in each layer were determined to obtain the best agreements for the temperature at the radius $r_0$ between measurements by DTS and calculations during heating and recovery. Theoretical calculations for soil temperature during recovery were conducted for constant and uniform heating [11]:

$$\Delta T_s\left(t>t^{\prime }\right)=\Delta T_s\left(t\right)-\Delta T_s\left(t-t^{\prime }\right),$$

(6)

where $t^{\prime }$ is the time at which heating stopped. The recovery period for analysis was defined as the time for the temperature change to drop to 10% or less its peak rate. $RMSE$ was used as a criterion for evaluating agreement between the calculated and measured temperatures [15]:

$$RMSE=\sqrt{\frac{1}{N}{{\displaystyle \sum }}_t{\left(\Delta T_{obs}-\Delta T_{cal}\right)}^2},$$

(7)

where $\Delta T_{obs}$ and $\Delta T_{cal}$ are the observed and calculated temperature changes from the initial temperature, respectively, N is the number (minute counts) of time steps for the history matching. $\Delta T_{obs}$ was obtained by averaging 42 × 2 DTS data at 0.5 m intervals in each 20.5 m thickness.

This study compared partial conductivities between ILS and MLS considering the uncertain effects of groundwater flow. First, the ILS theory in Equation (2) was used, resulting in partial estimates, including groundwater flow effects, i.e., the apparent thermal conductivity. To assess the effect of advection by groundwater flow, this study also applied another solution of the MLS theory instead of Equation (2) [44]:

$$\Delta T_s=\frac{q}{4\pi \lambda }{{\displaystyle \int }}_0^{\pi }exp\left(\frac{Ur}{2a}\cos \phi \right){{\displaystyle \int }}_0^{\frac{4at}{r^2}}\frac{1}{\eta }exp\left(-\frac{1}{\eta }-\frac{U^2{r}^2\eta }{16a^2}\right)d\eta d\phi ,$$

(8)

$$U=\frac{{\rho }_{\mathrm{w}}{C}_{\mathrm{w}}}{\rho C_{\mathrm{p}}}u$$

(9)

where $U$ is the modified velocity of groundwater [m/s], ${\rho }_{\mathrm{w}}{C}_{\mathrm{w}}$ is the heat capacity of groundwater (4.2 MJ/(m³·K)) and $\rho C_{\mathrm{p}}$ is that of the soil (3 MJ/(m³·K)), $u$ is the Darcy velocity (specific discharge) [m/s], and $\phi$ is flow direction (angle) [rad]. Integrating Equation (8) with respect to the flow angle $\phi$ results in the average temperature change at radius r_o. This study assumed a negligible effect of hydraulic dispersion under a rapid groundwater flow to simplify the calculation convergence, as applied in other studies [48]. Soil temperatures were calculated based on Equations (8) and (9), and then the vertical partial properties of $\lambda$, $q$, and $\log u$ were determined using the least $RMSE$ in each layer I, II, III, or IV at the eight BHEs. The Darcy velocity was converted to logarithmic form through parameter optimization, considering that the velocity was potentially variable by several orders of magnitude. The non-linear optimization algorithm “fmincon” in MATLAB was used to obtain the least $RMSE$, within the following range for each parameter: $1\le \lambda \ll 4,10\le q\ll 100,0.1\le \log u\ll 3$. Comparison of the parameters determined using ILS or MLS in each of the four layers aligning with the eight BHEs gives the site-scale variability with respect to the site’s geothermal and hydrogeologic properties.

4. Results

4.1. Conventional TRT Analysis

Table 1. Summary of configurations and analysis results for conventional Thermal response tests (TRTs).

BHE	U-Tube	Ls [m]	$\mathit{Q}$ [kW]	$\overline{\mathit{\lambda }}[\mathbf{m}·\mathbf{K}/\mathbf{W}]$	${\mathit{R}}_{\mathbf{b}}[\mathbf{m}·\mathbf{K}/\mathbf{W}]$
1	Cd		4.18	1.95	0.077
2	Cd+	82	4.42	2.17	0.089
3	C		4.01	4.93	0.114
4	C+	82	4.47	1.76	0.129
5	O		4.43	2.07	0.109
6	O+	10	4.01	1.87	0.110
7	O+	40	4.00	1.94	0.098
8	O+	82	4.05	1.83	0.106

summarizes the test configurations and analysis results, $\overline{\lambda }$ and $R_{\mathrm{b}}$, from conventional TRTs in the eight BHEs, and is modified from Table 1 in our previous paper [43] by the addition of the results for $\overline{\lambda }$. Conventional estimation of $\overline{\lambda }$ for the same length (82 m) ranged between 1.83 and 2.17, except for BHE No. 3; the average of the other seven BHEs was 1.94 W/(m·K). Considering a 10% error as ordinarily acceptable in conventional TRTs [17,49], the ground thermal conductivity in the seven BHEs was approximately equal to the average value, which itself was equal to reference values for dense gravel deposits [1,2,8]. However, at BHE No. 3, the estimate of $\overline{\lambda }$ was 4.93 W/(m·K), which was about 2.5 times greater than that at the other BHEs.

Figure 4 shows the changes of $\overline{T_f}$ measured by the PT100 sensors in the eight BHEs. $\overline{T_f}$ increased linearly in logarithmic time except at BHE No. 3. The intercepts of the fitted lines for $\overline{T_f}$ at the other seven BHEs depended on the heating rate in the two equipment sets, and the borehole thermal resistance related to the U-tube configuration. The plots for the seven BHEs showed slopes that were not clearly different from each other, indicating small variability of $\overline{\lambda }$. In contrast, $\overline{T_f}$ at BHE No. 3 also increased with time, but not linearly as at the other BHEs; it almost converged as the heating time progressed. The small increase in $\overline{T_f}$ was probably due to the effects of groundwater flow, as discussed later.

In-situ R_b values were compared and discussed in detail in the previous study [43]. This paper briefly summarizes the comparison. In Table 1, the in-situ R_b value of the oval/single(O) U-tube (No. 5) was about 10% lower than that of the circle/single(C) U-tube (No. 3), indicating the oval U-tube more effectively enhanced the thermal performance of the BHEs than the circular U-tube. However, the in-situ R_b value of oval/single tube (No. 5) was not lower than that of the more expensive circle/double (Cd) (No. 1), and thus, an oval-U tube is expected to be cost-effective relative to a double U-tube. The in-situ R_b value was not significantly different among the BHEs with and without spacers; Cd+ (No. 1) and Cd (No. 2), C (No. 3) and C+ (No. 4), O (No. 5) and O+ (Nos. 6, 7, and 8). The small effect of the spacer was due to the soil pressure in the deep zone; a part of each spacer was compressed with the buried pressure, and the spaces between the U-tubes were not maintained as expected. Future work will seek to develop the spacer to withstand soil pressure.

4.2. Vertical Profiles of DTS Measurements

Figure 5 shows vertical profiles for the initial soil temperature before heating $T_0$, and for the peak temperature after heating for 48 h, in each BHE. The initial temperature decreased with depth in all BHEs. However, No. 3 showed a different downward trend from the other seven BHEs. The seven BHEs had vertical decreases of T₀ with depth until about 50 m, below which it was almost stable, while No. 3 had a largely constant initial temperature between 10 and 50 m. This relative stability indicated active heat transfer, probably by groundwater flows.

The peak temperatures also showed different downward trends between BHE No. 3 and the other seven BHEs. The seven BHEs had peak temperatures above 20 °C at almost all depths. There was some small fluctuation with depth, as compared with the initial profiles, probably indicating that the uniformity of heat flux within the geologic structures at each depth. In contrast, the peak temperatures in BHE No. 3 were <20 °C at all depths, with the lowest value at shallow depth between 10 and 30 m. This means that groundwater flowed uniformly in the formation, and the flows were strong in the deep sections.

4.3. History Matching Analysis

Table 2. Summary of history matching analysis results.

BHE	Layer	ILS			MLS
		$\mathit{\lambda }$	$\mathit{q}$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{\lambda }$	$\mathit{q}$	$\mathbf{log}\mathit{u}$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$
1	I	2.15	49.39	0.28	2.11	49.26	1.41	0.28
	II	2.24	49.37	0.35	2.26	49.32	1.09	0.42
	III	2.34	49.37	0.39	2.41	49.27	1.19	0.46
	IV	2.51	49.40	0.45	1.59	34.91	1.60	0.35
2	I	2.00	49.38	0.20	1.24	35.51	1.38	0.12
	II	2.06	49.37	0.28	2.02	49.37	1.01	0.33
	III	2.27	49.39	0.37	1.32	34.58	1.60	0.22
	IV	2.38	49.39	0.43	2.25	49.03	1.66	0.39
3	I	4.72	49.27	0.51	2.39	42.71	2.24	0.16
	II	5.62	49.15	0.42	1.81	37.69	2.36	0.08
	III	3.75	49.41	0.45	2.13	39.68	2.01	0.15
	IV	3.66	49.38	0.33	2.09	37.06	1.87	0.18
4	I	2.02	63.86	0.76	1.73	64.17	1.74	0.32
	II	1.56	49.51	0.52	1.88	61.49	1.57	0.18
	III	1.65	49.51	0.54	1.48	49.53	1.63	0.19
	IV	1.71	49.48	0.79	1.44	49.58	1.75	0.22
5	I	2.16	49.38	0.28	2.22	49.23	1.34	0.26
	II	1.34	34.64	0.26	1.36	34.65	1.20	0.27
	III	2.12	49.36	0.37	1.39	34.65	1.30	0.29
	IV	2.26	49.41	0.46	1.49	35.03	1.66	0.20
6	I	2.07	49.40	0.15	1.98	49.35	1.31	0.12
	II	2.19	49.41	0.24	2.07	49.30	1.53	0.12
	III	2.18	49.39	0.27	2.06	49.30	1.56	0.15
	IV	2.21	49.40	0.30	2.09	49.41	1.60	0.15
7	I	2.26	49.40	0.24	2.25	49.15	1.48	0.16
	II	2.25	49.36	0.26	1.45	34.67	1.32	0.14
	III	2.34	49.36	0.29	2.28	48.88	0.18	0.32
	IV	2.46	49.39	0.40	1.62	35.01	1.71	0.13
8	I	1.86	49.40	0.20	1.79	49.36	0.61	0.17
	II	2.04	49.39	0.22	1.27	34.68	1.48	0.13
	III	2.09	49.41	0.48	2.00	49.41	1.72	0.18
	IV	2.20	49.41	0.35	2.17	49.31	1.64	0.18

summarizes the history matching analysis results: vertically partial estimates of ground thermal conductivity, heating exchange rate, and logarithmic Darcy velocity $\log u$, and the least $RMSE$. Figure 6 shows example results for BHEs Nos. 1 and 3 using both ILS and MLS solutions. For BHE No. 1, the estimates of ${\lambda }_{\mathrm{I}}$ were 2.15 and 2.11 W/(m·K) using ILS and MLS solutions, respectively. These values were approximately equal, but a comparison of Figure 2a and b shows that the fitting line more closely follows the MLS solution than the ILS solution. The MLS solution fitted more closely than the ILS solution most results in Table 2, as indicated by lower $RMSE$ values. The applicability of a parameter was dependent on the number of parameters, i.e., three in the MLS solution, but two in the ILS solution. However, except for BHE No. 3, the difference in $RMSE$ between the ILS and MLS solutions was small, within 0.15 K, indicating limited groundwater flow effect. However, in BHE No. 3, the MLS solution had much smaller $RMSE$ than the ILS solution. Figure 6 also indicates that the MLS solution provided much closer fitting than the ILS solution, especially for the stable temperatures after heating for 20–96 h. The clear applicability of the MLS solution means that the underground temperatures for history matching were affected significantly by groundwater flows and that the MLS solution was required to obtain a reasonable estimate of ground thermal conductivity in each layer.

Figure 7 shows box plots of partial $\lambda$ in the seven BHEs excluding No. 3, and individual plots for No. 3 using the ILS (a) and MLS (b) solutions. The boxes for the partial λ using the ILS solution ranged between 1.6 and 2.5 W/(m·K). The median values ranged between 2.1 and 2.3, with a slight increase in the deeper layers. The plots for estimates of partial λ in BHE No. 3 showed much larger medians in each layer than those in the other BHEs and those from conventional TRT analysis (Table 1). The plotted differences between BHE No. 3 and the others were typically larger in the shallower layers (I and II), corresponding to the depth sections with lower peak temperatures in Figure 5. The estimates of the partial ${\lambda }_{\mathrm{I}}, {\lambda }_{\mathrm{II}}, {\lambda }_{\mathrm{III}}, \mathrm{and} {\lambda }_{\mathrm{IV}}$ in BHE No. 3 were 4.7, 5.6, 3.8, and 3.7 W/(m·K), respectively, with an average of 4.45 W/(m·K), which was close (an error of about 11%) to the overall $\overline{\lambda }$ of 4.93 W/(m·K) in BHE No. 3.

The box plots of $\lambda$ for MLS solutions ranged between 1.4 and 2.2 W/(m·K) in layers I–IV, which were shifted to the left relative to the box plots for the ILS solutions. The median values for the four layers ranged from 1.6 to 2.0 W/(m·K). The narrow range was close to the overall average value for $\overline{\lambda }$ (1.94 W/(m·K)) at the seven BHEs. The estimates of the partial $\lambda$ in BHE No. 3 fell in or near the box plots for the other seven BHEs, meaning that the effect of groundwater advection was successfully extracted from the estimates of partial $\lambda$ in BHE No. 3.

Figure 7c shows box plots of $\log u$ for each layer in the seven BHEs and individual plots for No. 3. In seven BHEs except BHE No.3, the plots of $\log u$ ranged between 0.78 log m/y (6.0 m/y) and 1.75 log m/y (56 m/y) in layers I–IV. The values were larger than the criterion for advection effect (2 m/y) [37], and thus, the partial estimates of $\lambda$ using the ILS solution were higher than those using the MSL solution, due to the effect of groundwater flow. In BHE No. 3, with high apparent thermal conductivity in Figure 7a, the estimates of $\log u$ were 2.24 log m/y (174 m/y), 2.36 log m/y (228 m/y), 2.01 log m/y (101 m/y), and 1.87 log m/y (74 m/y) in layers I, II, III, and IV, respectively. History matching analysis based on the underground temperatures revealed rapid flows of groundwater in shallower layers I and II, only in the typical BHE among eight aligned BHEs.

5. Discussion

This study performed conventional TRT analysis for eight BHEs aligned on a line. The results for BHE No. 3 were typically higher than those for the others. History matching analysis was also performed based on the underground temperatures from DTSs. The MLS solution gave lower $RMSE$ than the ILS solution in most analysis cases, especially for BHE No. 3. The effectiveness of the MLS solution in the analysis indicates that groundwater flows affected the underground temperature during the DTRTs.

Using the MLS solution, median values of the partial $\lambda$ at seven BHEs and individual estimates for BHE No. 3 were obtained in a narrow range between 1.6 and 2.0 W/(m·K). The values were consistent with expectations for the gravel deposits that mainly comprised the formation. The median Darcy velocities in seven BHEs ranged between 6.0 and 56 m/y; those for BHE No. 3 ranged between 74 and 204 m/y. Considering the multiplications of the hydraulic conductivity of gravel deposits ranged between 10⁻³ and 10⁻⁰ m/s [50], and the topographic slope on the site was ranged between about 0.01 (1%) and 0.001 (0.1%), simple estimates for the Darcy velocity were ranged between 3.2 × 10¹ and 3.2 × 10⁵ m/y. The range of the simple estimates was reasonably agreed with the estimates by DTRTs, considering that the gravelly deposits were not clean and sorted by the inclusion of sand and clay deposits.

In conclusion, the ground thermal conductivity was affected by groundwater flows at each depth and location within this site, but the effect was mostly limited (i.e., at seven of the eight BHEs). On the other hand, rapid groundwater flows with Darcy velocity over 100 m/y were also possible, like at BHE No. 3. The rapid flows occurred in highly permeable zones within the gravel deposits. The gravel deposits around the foot of the mountains were not well-sorted, and porosity was distributed uniformly among them. The extremely high permeability was probably related to buried zones of open-void connectivity at shallow depths (less than 20 m) in BHE No. 3. The open-void connectivity was reported in the shallow zones of gravel deposits by sedimentary studies [51,52,53]. However, the permeable zone was limited to within a radius of a few meters from BHE No. 3.

Variability of ground thermal conductivity was observed across the study site. This field study provided a useful evidence of site-scale variability as an alarm for design and planning of GSHP systems at other sites on and around steep reliefs. High values of ground thermal conductivity were limited within this site, and similar could be expected at other sites. If a TRT detects a high ground thermal conductivity, the test BHE can be expected to obtain high thermal performance, and such a specific BHE of high ground thermal conductivity can be available in a residential or small office building of limited thermal loads. However, the high ground thermal conductivity might be not obtainable in other BHEs. In general, a high value in a steep relief should not be assumed representative for the design of multiple BHEs, considering the uncertainty in site-scale changes of ground thermal conductivity. Especially if groundwater flow simulation indicates rapid groundwater flows in some areas of steep relief, we should consider the uniformity of groundwater flows and take care not to overestimate advection effects.

Moreover, this study showed the effectiveness of the history matching analysis coupling with MLS to determine the vertical change in partial $\lambda$ affected by groundwater flows. The information of the vertical change was available to determine lengths of BHEs reasonably to reduce initial costs of the GSHP systems. However, the analysis required the measurements by DTSs, and the additional costs of DTSs might be no less than the cost reduction by optimal design of BHEs. In addition, the accuracy of the history matching is affected by the assumptions for theory of MLS. One assumption is the uniformity of the heat exchanger rate within each layer, and the layer span for the uniformity is uncertain. In this study, BHE was divided into four layers of the same length, in each of which the heat exchange rate was assumed to be constant during the testing. To address the problem, our previous study applied a pair of DTSs into the U-tubes to estimate the heat exchange rates in arbitrary depth sections, resulting in an optimal combination of partial conductivity and depth sections [27]. Besides, the analysis results included uncertainty, due to the measurement errors and the limitation of the assumptions. The uncertainty in Equation (4) for the ILS solution was derived as a sum of squared relative errors of measurements [11]. Therefore, the MLS solution’s uncertainty will be evaluated through sensitivity analysis of parameters in the history matching.

The MLS solution also neglected the change in thermal diffusivity, due to groundwater velocity. The hydraulic dispersion effect should be considered for more accuracy [54], although the uncertainty would increase by adding another parameter for optimization. In general, realistic results (in terms of minimized $RMSE$) might not result from using ILS of MLS. One suggestion to improve estimations using the ILS or MLS solutions is to repeat TRTs at different output powers $Q$, as the temperature changes are linear, due to heat conduction only, but non-linear when heat conduction occurs along with advection. Another approach is numerical modeling of detailed geologic and hydrogeologic conditions. These approaches will be demonstrated at the present site or at another with rapid groundwater flows in future work.

6. Conclusions

This study aimed to evaluate site-scale variability of ground thermal conductivity through TRTs in eight BHEs aligned on a line in steep relief with rapid and uniform groundwater flows, as a case study in Japan. The site was located on a terrace in the foothills of mountains in northern Japan. The formation consisted of monotonic gavel deposits. The BHEs were constructed with the same installation length (82 m) at intervals of 4 or 6.5 m. First, conventional analysis based on the temperature of circulating fluids during heating was performed to obtain overall values of ground thermal conductivity. Next, vertical partial values of ground thermal conductivity and Darcy velocity in four divided layers were determined through history matching analysis based on the underground temperatures measured by DTSs during heating and cooling. The results are summarized, as follows.

• The overall ground thermal conductivity was approximately equal at seven of the eight BHEs, with an average of 1.94 W/(m·K). BHE No. 3 had an overall value of 4.93 W/(m·K), about 2.5 times larger than the average of the other BHEs.
• History matching analysis showed closer-fitting temperature calculations using the MLS solution than the ILS solution. The improvement shown using the MLS solution was moderate in the seven BHEs, but clearer in No. 3. This study considered the extremely large value to be due to the effect of groundwater advection. The history matching coupling with the MLS solution was effective in evaluating the vertical change in the partial ground thermal conductivity affected by groundwater flows.
• The partial values of ground thermal conductivity using the MLS solutions were within a narrow range of 1.9–2.1 W/(m·K), which is reasonable for the site’s gravel deposits. The estimated Darcy velocities were 74–204 m/y at BHE No. 3, with the higher velocities seen in the shallow layers less than 40 m deep. The velocity at the other BHEs was one-to-two orders of magnitude smaller than that at BHE No. 3. These velocity values were reasonable for the ground elevation gradient and permeability of the gravel deposits.
• The high ground thermal conductivity at BHE No. 3 was limited to within a radius of a few meters. Such a spatial variability of ground thermal conductivity is considered in any other sites on and around high reliefs. This means that the effects of advection on the thermal performances of BHEs in the steep ground (even if the high value was actually obtained by TRT) should be carefully considered when designing and planning GSHP systems.