A New Estimate of Sand and Grout Thermal Properties in the Sandbox Experiment for Accurate Validations of Borehole Simulation Codes

Abstract

Ground-coupled heat pumps usually employ fields of borehole heat exchangers (BHEs), which must be designed by suitable models. In order to validate a BHE model, it is advisable to compare the computation results with experimental data. A well-known data set was provided by Beier et al. (Geothermics 2011, 40) through a laboratory model usually called “sandbox”. Several authors proposed estimates of the thermal properties of the sandbox grout and sand. In this paper, we present a new estimate of those properties, obtained by means of 2D finite-element simulations that consider all the details of the experimental setup, including the thin aluminum pipe at the BHE boundary. Our results show that the measured temperatures in the fluid and in the sand can be reproduced very accurately by considering thermal conductivities 0.863 W/(mK) for the grout and 3.22 W/(mK) for the sand, volumetric heat capacities 4.6 MJ/(m³K) for the grout and 3.07 MJ/(m³K) for the sand, and a slightly enhanced heat capacity of the water contained in the BHE. The 2D simulations are validated by comparison with an analytical solution and by 3D simulations.

1. Introduction

Applications of ground-source heat pumps (GSHPs) have rapidly spread during the last years since they represent an efficient technology for building heating and cooling [1]. GSHPs include ground-coupled heat pumps (GCHPs) [2], ground-water heat pumps (GWHPs) [3,4,5,6] and surface-water heat pumps (SWHPs) [7]. GCHPs usually make use of vertical ground heat exchangers, also named borehole heat exchangers (BHEs), which are composed of either a single or a double U-tube connected to the ground by a suitable grouting material. Analytical or numerical BHE models are needed for the design of BHE fields, for the analysis of thermal response tests (TRTs) [8], and for the simulation of GCHP systems [9].

The simplest BHE models are the infinite line-source (ILS) model and the finite line-source (FLS) model, where the BHE is sketched either as a line placed in an infinite solid medium or as a line segment placed in a semi-infinite solid medium. Contrary to the ILS model, the FLS one considers the heat transfer in the vertical direction and can be also employed for long-term simulations of a BHE. Analytical expressions of the temperature field produced by a FLS that releases a uniform and constant heat rate per unit length have been developed by Claesson and Eskilson [10] and by Zeng et al. [11], in the case of a line-source with the top at the ground surface. Simplified expressions of the temperature field averaged along the BHE length have been determined by Lamarche and Beauchamp [12] and by Bandos et al. [13]. A solution for the more general case in which the BHE top is buried under the ground surface has been presented by Claesson and Javed [14]. The line-source models yield accurate values of the mean temperature of the BHE surface after some hours from the beginning of the heat transfer process, but are not precise during the first few hours. Moreover, they require the use of the BHE thermal resistance, R_b, to determine the mean fluid temperature, and R_b is a well-defined property of the BHE geometry and materials only when the heat transfer in the BHE is quasi-stationary, i.e., after some hours of operation. Therefore, more accurate BHE models have been developed to determine the time evolution of either the mean fluid temperature or the fluid temperature at the BHE outlet in the short term.

Capacity-resistance BHE models, which are also accurate in the short term, were developed by De Carli et al. [15], Zarrella et al. [16], Bauer et al. [17], Pasquier and Marcotte [18,19], Ruiz-Calvo et al. [20], and Mineai and Maerefat [21]. Li and Lai [22,23] determined the analytical solution for a 2D short-term BHE model in which the pipes are represented by infinite line sources. Zhang et al. [24] proposed a quasi-3D line source model that employs transient BHE thermal resistance and holds for the entire time scale.

Beier and Smith [25] found an analytical solution for a cylindrical BHE model in which the real BHE is replaced by a grout annulus with external radius equal to the BHE radius and the same thermal resistance as the BHE. Xu and Spitler [26] developed a numerical BHE model where all the materials of the BHE, including the fluid, are represented by concentric cylindrical layers. Man et al. [27] proposed a simple analytical BHE model, where the BHE materials are replaced by the ground and the fluid is represented by a heat generating cylindrical surface.

Bandyopadhyay et al. [28] proposed an analytical BHE model where the fluid is represented by a heat-generating solid cylinder with high thermal conductivity, surrounded by a grout layer. The authors only determined the analytical solution in the Laplace transformed domain, so that a numerical inversion was required. A complete analytical solution for a BHE model similar to that proposed in Reference [28], with heat generation replaced by an injected heat flux, was given by Javed and Claesson [29].

Beier [30] determined the analytical solution, in the Laplace transformed domain, of a BHE model where the U-tube is represented by two half tubes surrounded by a cylindrical grout layer. Lamarche [31] found the analytical solution for a BHE model where the fluid is represented by a heat-generating solid cylinder surrounded by cylindrical layers that represent the tubes and the grout. Naldi and Zanchini [32] proposed a numerical model where the BHE is represented by a one-material equivalent cylinder (OMEC) with an internal heat-generating surface that represents the fluid. The authors evaluated the time evolution of the mean fluid temperature for a square field composed of four BHEs fed in parallel [33] by employing the OMEC model and an iterative method to impose a uniform temperature for the BHE fluid [34].

In order to check the accuracy of a BHE model, it is necessary to compare the computation results with those obtained by another model, or better yet, with experimental data. Several field tests on real BHEs have been performed, but these tests generally lack independent measurements of soil properties [35] or temperature measurements in the soil. Additionally, the position of the U-tube in the borehole is difficult to control and detect in the field [35]. In order to perform a rigorous experimental validation, reference data sets obtained by laboratory models are needed, with measured time evolutions not only of the inlet and the outlet fluid temperature, but also of the temperature of the BHE surface and of the surrounding soil.

An excellent and well-known reference data set has been provided by Beier et al. [35], who performed a thermal response test on a large laboratory BHE model and measured the time evolution of the inlet and outlet fluid temperature, of the BHE surface temperature at four axial locations, and of the surrounding sand at four axial locations and four radial distances from the BHE surface. The thermal conductivity of the grout was measured in a portion of the grout mixture, placed in a container, by means of a non-steady-state thermal probe. The thermal conductivity of the sand was measured at four locations along the BHE, by inserting a non-steady thermal probe in the sandbox. The results were 0.73 W/(mK) for the grout and 2.82 W/(mK) for the sand. The thermal conductivity and the volumetric heat capacity of the sand were estimated later by Beier [30] as 2.94 W/(mK) and 3.2 × 10⁶ J/(m³K), by analyzing the heat conduction in the sand.

Several authors employed the sandbox experiment by Beier et al. [35] to validate their BHE heat transfer models. Some of them also provided new estimates for the properties of grout and ground.

Pasquier and Marcotte [19] validated their quasi-3D thermal resistance capacity model (TRCM) for the simulation of BHEs by comparison with the data set provided by Beier et al. [35]. For the validation, they employed a thermal conductivity of 2.88 W/(mK) and a volumetric heat capacity of 2.55 MJ/(m³K) for the sand, and a thermal conductivity of 0.90 W/(mK) for the grout. The latter is much higher than that reported in Reference [35], and was obtained by imposing that the BHE thermal resistance was equal to that measured in Reference [35]. The authors obtained a very good agreement with the experimental evolution of the inlet and of the outlet fluid temperature.

Minaei and Maerefat [21] validated their 3D TRCM model by comparison with the data set provided by Beier et al. [35]. They considered the same geometric parameters and the same thermal properties reported by Beier [30], except for grout and ground thermal conductivities. For these quantities, the authors considered the values estimated by Pasquier and Marcotte [19].

Boockmeyer and Bauer [36] performed accurate 3D numerical simulations of the laboratory BHE model in Reference [35], and found that the input parameters given in Reference [35] yield inlet and outlet fluid temperatures higher than the measured values. Therefore, they performed an estimate of the grout and sand thermal properties and of the thermal power in the sandbox experiment [35], through their 3D simulation code and the parameter estimation software PEST [37]. They found the following results: thermal conductivity and volumetric heat capacity of the grout 0.80 W/(mK) and 3.80 MJ/(m³K), thermal conductivity and volumetric heat capacity of the ground 2.80 W/(mK) and 2.92 MJ/(m³K), thermal power 1013 W. The error in the measurement of the electric power estimated in Reference [36], + 4.1% is much higher than that declared in Reference [35], ±1%.

Zhang et al. [24] employed the reference data set of Reference [35] to validate their quasi-3D line source model of BHEs, which employs a transient BHE thermal resistance. They found that the grout thermal conductivity given in Reference [35] yields both a BHE thermal resistance higher than that measured directly by Beier et al. [35] and an overestimation of the inlet and outlet fluid temperatures, while the use of a modified transient BHE thermal resistance based on the value of steady-state thermal resistance measured directly in Reference [35] yields the correct time evolution of the fluid temperatures.

Li and Lai [23] validated their infinite composite-medium line-source (ICMLS) analytical model for the short-term response of U-tube BHEs by comparison with the data set in Reference [35]. The model [22,23] replaces the legs of a U-tube by two infinite line sources, each placed in the center of a leg, and considers the different thermal properties of grout and ground. The authors found that by applying to the model the input data ofReference [35], one finds an overestimation of the inlet and of the outlet fluid temperatures with respect to the measured values. An overestimation of the outlet fluid temperature with respect to the measured values of Reference [35] was also found by Ma et al. [38], by employing an analytical BHE model based on that developed in References [22,23], completed for applications in the long term.

Zhu et al. [39] and Lei et al. [40] validated their 3D numerical BHE models by comparing the computed time evolution of the outlet fluid temperature, of the temperature of the BHE surface and of the temperature of the sand at different radial distances with the experimental data in Reference [35]. The authors employed the inlet fluid temperature as an input, and obtained results in fair agreement with the experiment. Similarly, Tombarević and Vušanović [41] and Rohde et al. [42] successfully validated their quasi 3D BHE models through a comparison to the reference data set in Reference [35] by assuming the inlet fluid temperature as an input. Likewise, Brussieux and Bernier [43] successfully validated a hybrid numerical/analytical one-dimensional BHE model, based on the equivalent geometry proposed by Xu and Spitler [26], by means of the reference data set provided by Beier et al. [35]. The authors took the system parameters, geometry and thermal conductivities from Reference [35], and the heat capacities from Minaei and Maerefat [21], and used the experimental values of flow rate and inlet temperature as inputs. Indeed, assuming the inlet temperature as an input dramatically reduces the effects of an incorrect estimation of the grout thermal conductivity on the simulation results. In fact, for a given thermal power supplied to the BHE, a lower (higher) thermal conductivity of the grout yields an increase (a decrease) in the fluid inlet temperature.

Cherati and Ghasemi-Fare [44] used the data set of Reference [35] to validate their analytical model of heat and moisture transfer in an unsaturated porous medium surrounding a heat source. They assumed the temperature rise at the borehole surface as an input and did not analyze the whole heat transfer problem in the BHE.

Li et al. [45] analyzed the heat transfer processes in the BHE system considered in Reference [35] by adopting two computational fluid dynamics models. In the first model, the authors assumed the measured inlet fluid temperature as an input condition; in the second, they assumed the measured power supplied to the BHE as an input. While the first model yielded a time evolution of the outlet fluid temperature in fair agreement with the experimental results, the second yielded an overestimation of the outlet fluid temperature up to 2 °C.

Wang et al. [46,47] validated by the data set of Reference [35] two infinite composite-medium cylindrical source (ICMCS) BHE models, based on the equivalent geometry proposed by Lamarche [31]. In these models, the BHE is composed of an inner cylinder, representing the fluid, surrounded by an annulus representing the equivalent pipe, surrounded in turn by an annulus representing the grout and having an external radius equal to the BHE radius. The external radius of the equivalent pipe is determined by imposing that the grout annulus has the same thermal resistance as that of the grout in the real BHE. The authors found good agreement with the experimental values of the inlet and outlet fluid temperatures. However, they assumed that the BHE thermal resistance is 0.158 mK/W, which corresponds to a grout thermal conductivity much higher than that measured in Reference [35].

Zhang et al. [48] applied to the data set in Reference [35] their method for the estimation of grout and soil thermal properties in TRTs, based on the ICMLS model. They assumed, as input parameters, the geometric data, the thermal conductivity of the U-pipe, the undisturbed ground temperature, the heat flux and the volume flow rate, and estimated the thermal properties of grout and ground. By considering the shank spacing measured in Reference [35], they found a thermal conductivity of the grout equal to 0.89 W/(mK), i.e., 21.9% higher than that measured by Beier et al. [35], and a thermal conductivity of the sand equal to 3.244 W/(mK), i.e., 15.04% higher. In order to reduce the discrepancy on the grout conductivity, they conjectured that the real shank spacing was higher than the measured one, namely 6.88 cm, instead of 5.3 cm. Zhang et al. [49] again analyzed the data set of Reference [35] by a modified version of the ICMLS model, under the assumption that the shank spacing and the thermal conductivities of grout and sand are those reported in Reference [35], but the pipes are not placed in correspondence to the BHE axis and touch the BHE surface. By this assumption, they reproduced the time evolution of the inlet and outlet fluid temperature with a maximum error of 0.8 °C, starting from 1 h after the TRT start.

Zhang et al. [50], by applying a parameter estimation algorithm based on the ICMLS model to the data set in Reference [35], found a value of the thermal conductivity of the grout slightly more than 16% higher than that measured in Reference [35].

To summarize, some authors [23,24,38,45] pointed out that the parameters reported in Reference [35] yielded a relevant overestimation of the time evolution of the inlet and the outlet fluid temperatures. Others [19,21,46,47,50] estimated higher values of the grout thermal conductivity or employed values of the BHE thermal resistance corresponding to higher values of the grout thermal conductivity, to obtain the correct time evolution of the fluid temperatures. Others [36] estimated both a higher grout conductivity and a lower thermal power released to the BHE. Some authors [39,40,41,42,43] found good agreement between the simulated and the measured time evolution of the outlet fluid temperature, with the parameters reported in Reference [35], by assuming the inlet fluid temperature as an input. This kind of input strongly reduces the effects of an incorrect estimation of the grout thermal conductivity. Finally, some authors [48,49], as an alternative to estimating a grout thermal conductivity much higher than that measured in Reference [35], conjectured that either the shank spacing was higher than that declared in Reference [35], or the pipes were not lying in the same plane as the BHE axis and touched the BHE surface.

In this paper, by means of 2D finite-element simulations, we provide new best estimates of the grout and of the ground thermal properties by a method that allows for separate estimations of these properties. In our analysis, we consider the presence of a thin aluminum pipe at the BHE boundary and the difference between the local temperatures of the BHE surface and of the sand measured in Reference [35], in correspondence with the inlet pipe, and the mean temperatures at the same radial distances.

2. The Sandbox Experiment

Beier et al. [35] created a large laboratory BHE model, composed of a U-tube in high density polyethylene surrounded by a bentonite grout. The U-tube had a pipe outer diameter of 3.34 cm, a pipe inner diameter of 2.733 cm, a shank spacing of 5.3 cm, and a length of 18.32 m. The grout was contained within an aluminum pipe with an external diameter of 13 cm and a thickness of 0.2 cm. The U-tube was centered in the borehole with spacers so that the same separation distance existed between the outer surfaces of the two pipes and the outer surface of each pipe and the borehole wall [35]. The BHE was inserted in a square wooden box with sides of 1.83 m, filled with wet sand. The ends of the wooden box were insulated. A stream of conditioned air was circulated over the lateral surface of the box in order to keep the temperature of this surface constant. A TRT with constant power supplied to the BHE and a duration of 3106 min was performed. The inlet and outlet water temperatures, the temperatures of the outer surface of the aluminum pipe, and the temperatures of the sand at distances of 24, 44, 65 and 85 cm from that surface (i.e., at radial distances 30.5, 50.5, 71.5 and 91.5 cm from the BHE axis) were measured at time intervals of 1 min by thermistors with an uncertainty of ±0.03 °C. A scheme of the experimental setup is reported in Figure 1, and some images can be found in Reference [35]. Figure 2 illustrates part of a cross section of the apparatus, showing the radial positions of the temperature sensors S1, S2, S3, S4, and S5. Each temperature sensor represented in Figure 2 corresponds to four thermistors of the experimental apparatus, placed at different axial positions, as shown in Figure 1. For each temperature sensor, we have considered the average of the temperatures measured by the corresponding thermistors.

The supplied electric power was measured at every minute with an uncertainty of ±1%; the mean value during the whole TRT was 1051.6 W. The flow rate, measured at every minute by a flow meter, had a mean value 0.197 L/s. The thermal conductivity of polyethylene pipes was declared to be k_p = 0.39 W/(mK). The thermal conductivity of the grout was measured by a non-steady-state thermal probe, inserted in a portion of the grout mixture placed in a container and allowed to cure. The measurement result was k_g = 0.73 W/(mK). The thermal conductivity of the sand was measured in a similar way, by a probe inserted directly in the sand through ports that could be opened and closed. The measurement result was 2.82 W/(mK). The volumetric heat capacity of the grout, (ρc)_g, and of the sand, (ρc)_s, were estimated later by Beier [30] as 3.8 MJ/(m³K) and 3.2 MJ/(m³K), respectively. All the measurement results of the TRT were kindly provided to us by Prof. Beier, who also suggested considering a thermal conductivity equal to 177 W/(mK), a density equal to 2770 kg/m³ and a specific heat capacity equal to 875 J/(kgK) for the aluminum pipe.

3. Estimation Method and Results

The new estimate of the thermal properties of the grout and of the sand performed in this paper is based on the following assumptions.

• The measurements of temperature and of the thermal power released to the BHE are precise. This assumption is based on the accuracies declared in Reference [35] and on the nearly perfect agreement between the value of thermal power released to the BHE determined by measurements of voltage and current and that determined by the energy balance equation in the fluid, in quasi-stationary working conditions. The discrepancy between the mean values of these quantities, in the time period between 1200 min and the end of the TRT, is lower than 0.5%.
• The most probable positions of the U-tube pipes are those declared by the authors [35], with axes of the U-tube pipes and axis of the aluminum tube laying in the same horizontal plane, and a shank spacing 5.3 cm. Indeed, the BHE was only 18.32 m long, and spacings were inserted to ensure that the distance between the pipes and the distance between the outer surface of the pipes and the BHE wall were uniform [35].
• Since the U-tube is very short, the bulk fluid temperature along the flow can be considered as a linear function of the distance from the inlet section. Therefore, the mean fluid temperature in the whole U-tube, T_fm, that in the inlet tube, T₁, and that in the outlet tube, T₂, can be evaluated as

  $$T_{fm}=\frac{T_{in}+T_{out}}{2}\hspace{0.33em},\hspace{1em}{T}_1=\frac{T_{in}+T_{fm}}{2}\hspace{0.33em},\hspace{1em}{T}_2=\frac{T_{out}+T_{fm}}{2}$$

  (1)

  where T_in is the inlet temperature and T_out is the outlet temperature.
• There may be a significant discrepancy between the value of the grout thermal conductivity, k_g, measured in Reference [35] and the real value, because the measurement was performed in a sample of the grout placed in a container and allowed to cure. The BHE grout, initially prepared in the same way, may have absorbed water from the surrounding wet sand, and an increase in water content may have produced an increase in k_g with respect to the measured value.

The method used in this paper to estimate the thermal properties of the grout and of the sand employed in the sandbox experiment is based on 2D finite-element simulations of a cross section of the apparatus, implemented in COMSOL Multiphysics. The software employs a finite-element method for the solution of partial differential equations, where the time discretization is based on the Backward Differentiation Formula (BDF) solver and the space discretization is performed by quadratic Lagrange polynomials. We selected unstructured triangular meshes, which more accurately subdivide very thin elements of the domain, such as the aluminum tube. Our estimation method is composed of eight parts, as described below.

3.1. Preliminary Estimate of k_g, with k_s = 2.82 W/(mK)

The experimental data allow for a direct measurement of the ratio

$$R_{b1}=\frac{T_{fm}-T_{s1}}{q_b}$$

(2)

where T_s1 is the temperature measured by the sensor S1, i.e., the temperature of the BHE surface averaged along the BHE length and measured at the aluminum external surface in correspondence with the center of the inlet tube, and q_b is the thermal power per unit length released to the BHE. Note that R_b1 does not correspond exactly to the BHE thermal resistance, R_b, defined as

$$R_b=\frac{T_{fm}-T_b}{q_b}$$

(3)

where T_b is the mean temperature of the BHE surface. In fact, the BHE surface temperature in correspondence with the tubes is higher than that at a higher distance from the tubes, and the mean temperature of the fluid in the inlet tube, T₁, is higher than that of the fluid in the outlet tube, T₂.

We determined our best estimate of k_g by imposing that the value of R_b1 obtained by a stationary 2D simulation of a BHE cross section equals that obtained by employing the experimental data. In the evaluation of R_b1 through the experimental data, we assumed q_b equal to the mean value during the whole TRT, namely q_b = 57.401 W/m (1051.6 W and 18.32 m), and T_fm–T_s1 equal to the mean value in the time interval from 1500 min to the end of the TRT, namely 9.919 °C. The result is R_b1 = 0.1728 mK/W. In the evaluation of R_b1 through the simulation, we considered a BHE cross section, including the surrounding sand. We considered the thermal conductivity of the sand given by Beier et al. [35], k_s = 2.82 W/(mK), and imposed the temperature 0 °C at the external boundary. We imposed different mean temperatures of the fluid in the inlet and in the outlet pipe, T₁ and T₂, such that the difference in T₁–T₂ was equal to 0.5 (T_in–T_out) and the value of q_b resulting from the simulation was nearly equal to the experimental value. We considered the mean value of T_in–T_out during the time interval from 1500 min to the end of the TRT, which yielded T₁–T₂ = 0.6316 °C, and the convection coefficient h = 1948.8 W/(m²K).

The convection coefficient was calculated through the Churchill correlation with uniform wall heat flux [51], by considering the thermal properties of water at 30 °C, taken from the NIST Chemistry WebBook [52], namely density 995.65 kg/m³, dynamic viscosity 0.79735 mPa s, thermal conductivity 0.61550 W/(mK), and specific heat capacity at constant pressure 4179.8 J/(kg K)). The values of the time-averaged volume flow rate, of the Reynolds, Prandtl and Nusselt numbers, and of the convection coefficient are reported in

Table 1. Values of volume flow rate, of Re, Pr, and Nu numbers, and of the convection coefficient.

Flow Rate [L/s]	Re	Pr	Nu	h [W/(m2K)]
0.197	11,460	5.4147	86.533	1948.8

.

Stationary 2D simulations of a cross section were performed through COMSOL Multiphysics by selecting relative accuracy 10⁻⁶ and by employing an unstructured mesh with 50,872 triangular elements. A particular of the mesh is illustrated in Figure 3. The simulation with k_g = 0.862 W/(mK), performed with T₁ = 19.22 °C, T₂ = 18.5884 °C, and as a consequence, T_fm = 18.9042 °C, yielded T_s1 = 8.98365 °C, q_b = 57.396 W/m, and R_b1 = 0.1728 mK/W, in perfect agreement with the experimental value. The simulation was repeated with a mesh obtained by a regular refinement of the previous one and having 203,488 triangular elements. The results coincided with those reported above and revealed that refinements of the adopted mesh are useless. To check the mesh independence of the results, the simulation was repeated with five different meshes, obtained by applying a sequence of regular refinements to a mesh with 5214 triangular elements, denoted by Mesh 1. The results obtained for R_b1 are reported in Figure 4, and show that the value of R_b1 obtained with the adopted mesh is extremely close to those obtained by Meshes 4 and 5. Therefore, the mesh independence is proved and the result of the preliminary estimate of the thermal conductivity of the grout is k_g = 0.862 W/(mK).

3.2. Preliminary Estimate of the Thermal Diffusivity of the Sand, α_s

To estimate the thermal diffusivity of the soil, α_s, we performed 2D simulations of the transient heat conduction in a cross section of the sand. In order to employ a uniform initial condition on the whole computational domain, in all our transient simulations we replaced the experimental values of the temperature, T(t), with those of the temperature rise with respect to the initial value, θ(t), defined as

$$\theta \left(t\right)=T\left(t\right)-T\left(t=0\right)$$

(4)

and we evaluated θ(t) in our simulations.

In the preliminary estimate, we considered the temperature rise of the BHE surface as uniform and equal to the temperature rise measured by S1, θ_s1(t). We imposed θ_s1(t) as a boundary condition at the external surface of the aluminum pipe, and considered the external surface of the sand, namely the square with a side measuring 183 cm and a center in the axis of the BHE as adiabatic. In fact, the mean temperature of that surface did not remain constant during the TRT, as planned in the experiment, and the adiabatic condition reproduced the increasing trend of the mean temperature of that surface rather well. We determined α_s by minimizing the root mean square deviation (RMSD) between the experimental values of the temperature rises measured by S2, S3, and S4, i.e., θ_s2(t), θ_s3(t), θ_s4(t), and those obtained by the simulations. We performed transient 2D simulations by COMSOL Multiphysics, with the solver PARDISO, absolute tolerance 10⁻⁴, and an unstructured mesh composed of 28,320 triangular elements, which is illustrated in Figure 5.

The lowest RMSD between the values of θ_s2(t), θ_s3(t), θ_s4(t) resulting from simulations and the experimental ones was obtained with α_s = 1.00 × 10⁻⁶ m²/s, and was equal to 0.0260 °C. The simulation with this value of α_s was repeated with a mesh obtained by a regular refinement of the previous one, having 113,280 triangular elements. The refined mesh yielded values of θ_s2(t), θ_s3(t), θ_s4(t) nearly identical to those obtained with the previous mesh, and an RMSD equal to 0.0264 °C Therefore, refinements of the adopted mesh are unnecessary. To check the mesh independence of the results, the simulation was repeated with four different meshes, obtained by applying a sequence of regular refinements to a mesh with 1130 triangular elements, denoted as Mesh 1. The values of θ_s2, θ_s3, θ_s4 obtained with Meshes 1, 2, 3, 4 and with the adopted mesh, for two selected time instants, are reported in

Table 2. Preliminary estimate of α_s: values of θ_s2, θ_s3, and θ_s4 at selected instants of time, obtained by different meshes.

Mesh	Elements	t = 1500 min			t = 3106 min
		θs2 °C	θs3 °C	θs4 °C	θs2 °C	θs3 °C	θs4 °C
1	1130	1.4529	0.5096	0.1592	2.4186	1.2300	0.6676
2	4520	1.4542	0.5105	0.1592	2.4198	1.2308	0.6676
3	18,080	1.4541	0.5104	0.1592	2.4198	1.2308	0.6676
4	72,320	1.4541	0.5104	0.1592	2.4198	1.2307	0.6676
Adopted	28,320	1.4541	0.5105	0.1593	2.4198	1.2308	0.6676

. The table shows that the values obtained by the adopted mesh coincide with those obtained by Meshes 3 and 4, except, in some cases, for 1 unit in the fourth decimal place. The convergence of Meshes 1, 2, 3, and 4 to the values yielded by the adopted mesh, for θ_s2, is illustrated in Figure 6. Therefore, the mesh independence is proved and the result of the preliminary estimate of the thermal diffusivity of the sand is α_s = 1.00 × 10⁻⁶ m²/s.

3.3. Preliminary Estimate of k_s

In order to estimate k_s, we performed 2D transient simulations of a cross section of the whole system BHE + soil. We employed the value of (ρc)_g given by Beier [30], namely 3.8 MJ/(m³K), and the values of k_g and of α_s determined in the previous simulations, k_g = 0.862 W/(mK) and α_s = 1.00 × 10⁻⁶ m²/s. Since, in a 2D simulation of a BHE cross section, it is impossible to represent the fluid flow along the pipes, we considered water as a solid with thermal conductivity 1000 W/(mK) and a uniform heat generation in each tube. Indeed, replacing the fluid by a heat generating solid is a well-established technique in BHE simulation models [28,31,33]. The mean value of the heat generation was set as equal to the ratio between the thermal power released by the BHE and the volume of water in the BHE, namely q_g = 48,924 W/(m³). The heat generation in the inlet tube was set as equal to 1.0364 q_g, and that in the outlet tube was set as equal to 0.9636 q_g, to obtain a difference θ₁–θ₂ at the final instant very close to the experimental value, namely 0.508 °C, given by one half of θ_in–θ_out. A time delay of 1 min in the heat generation was imposed in the outlet tube, in order to take into account the flow time of water in the inlet tube. The time delay was implemented by employing an increasing function with continuous first and second derivatives, between 30 s and 90 s, centered at t = 60 s. To take into account the thermal inertia of portions of the circuit external to the BHE, the total volume of water in the system was assumed to be equal to that contained in the BHE tubes plus 8 L of external water. Thus, the density of the solid representing water was set as equal to that of water multiplied by 1.3722. The convective thermal resistance between the fluid and the inner surface of the pipes was inserted in the model as an increased thermal resistance of the pipes. Therefore, the thermal conductivity of the pipes, k_p = 0.39 W/(mK), was replaced by an effective thermal conductivity, which embodies the convective thermal resistance: the obtained value is k_peff = 0.36346 W/(mK). The external surface of the computational domain was considered as adiabatic. The value of k_s was estimated by minimizing the mean square deviation between the experimental values of θ₁(t), θ₂(t), θ_s1(t), θ_s2(t), θ_s3(t), θ_s4(t) and those obtained through the simulations. The simulations were performed by employing the same solver and absolute accuracy, as in Section 3.3, and an unstructured mesh with 52,152 triangular elements. A particular of the mesh is illustrated in Figure 7.

The lowest RMSD between the values of θ₁(t), θ₂(t), θ_s1(t), θ_s2(t), θ_s3(t), θ_s4(t) resulting from simulations and the experimental ones were obtained with k_s = 3.18 W/(mK), and was equal to 0.1308 °C. A simulation with a mesh obtained by a regular refinement of the previous one, with 208,608 triangular elements, yielded nearly identical results, and the same RMSD. Therefore, refinements of the adopted mesh are unnecessary. The mesh independence was verified by repeating the simulation with three different meshes, with 5472, 21,888 and 87,552 triangular elements. The values of θ₁, θ₂, θ_s1 obtained by these meshes, for t = 1500 min and t = 3106 min, are reported in

Table 3. Preliminary estimate of k_s: values of θ₁, θ₂, and θ_s1 at selected instants of time, obtained by different meshes, with k_g = 0.862 W/(mK), (ρc)_g = 3.8 MJ/(m³K), α_s = 1.00 × 10⁻⁶ m²/s, and k_s = 3.18 W/(mK).

Mesh	Elements	t = 1500 min			t = 3106 min
		θ1 °C	θ2 °C	θs1 °C	θ1 °C	θ2 °C	θs1 °C
1	5472	15.7905	15.2773	5.6623	16.9054	16.3921	6.7453
2	21,888	15.7912	15.2780	5.6626	16.9061	16.3928	6.7456
3	87,552	15.7912	15.2780	5.6626	16.9062	16.3929	6.7456
Adopted	52,152	15.7912	15.2780	5.6626	16.9062	16.3929	6.7456

. The table shows that the values obtained by the adopted mesh coincide with those obtained by Meshes 2 and 3, so that the mesh independence is proved and the result is k_s = 3.18 W/(mK).

3.4. Estimate of k_g, with k_s = 3.18 W/(mK)

Since the thermal resistance of a BHE depends, although slightly, on the thermal conductivity of the surrounding soil, we repeated the estimate of k_g described in Section 3.1 with the same method and mesh, and the value of k_s determined in Section 3.3, namely 3.18 W/(mK). The simulation with k_g = 0.863 W/(mK), performed with T₁ = 18.22 °C, T₂ = 17.5884 °C, and as a consequence, T_fm = 17.9042 °C, yielded T_s1 = 7.9846 °C, q_b = 57.417 W/m, and R_b1 = 0.1728 mK/W, in perfect agreement with the experimental value. We repeated the simulation with a regularly refined mesh with 203,488 triangular elements, and we obtained the same results. Due to the very high similarity with the conditions considered in Section 3.1, further checks of mesh independence did not seem necessary. Therefore, our estimate of the thermal conductivity of the grout is k_g = 0.863 W/(mK).

3.5. Estimate of the Time Evolution of the Temperature Distribution along the BHE Surface

The preliminary value of α_s was determined by considering the temperature rise of the BHE surface as uniform and equal to θ_s1(t). In order to improve the estimate, we imposed, along the BHE surface, the time-dependent distribution of θ (t) obtained by repeating the final simulation of Section 3.3 with the thermal conductivity of the grout k_g = 0.863 mK/W, determined in Section 3.4. To determine the expression of θ (t), we computed the temperature-rise distribution of the upper part of the BHE surface at the final time instant, as a function of the arc length, starting from the point opposite the sensor that measured θ _s1. We performed a fourth-order polynomial best fit of that distribution, as a function of the axial coordinate, x. The result was

$$\theta \left(x\right)={\theta }_{s1}+3305.3{\left(x-0.065\right)}^4+ 873.12{\left(x-0.065\right)}^3+ 130.96{\left(x-0.065\right)}^2+ 9.9823\left(x-0.065\right)$$

(5)

Then, we computed the time evolution of the ratio, A(t), between the difference θ _s1–θ _s0 at the final instant and that at the instant t, where θ _s0 is the temperature rise of the highest point of the BHE surface, namely that in the upper part for x = 0. Finally, we expressed the time-dependent temperature-rise distribution on the BHE surface by the equation:

$$\theta \left(x,t\right)={\theta }_{s1}\left(t\right)+A\left(t\right)\left[\begin{array}{l}3305.3{\left(x-0.065\right)}^4+ 873.12{\left(x-0.065\right)}^3+ 130.96{\left(x-0.065\right)}^2\\ + 9.9823\left(x-0.065\right)\end{array}\right]$$

(6)

Several numerical checks revealed that Equation (6) becomes very accurate when the difference θ _s1–θ _s0 becomes relevant, so that it gives an excellent representation of the time-dependent temperature distribution along the BHE surface for our purposes. In Figure 8, the temperature rise for the final instant, θ _fin, and that for t = 100 min, θ ₁₀₀, are plotted versus the angle, ψ, along the upper part of the BHE surface, starting from the point opposite S1. The figure reports both the values obtained by a 2D finite-element simulation and those determined by Equation (6), denoted by θ _{fin_approx} and θ _{100_approx}. The 2D simulation was performed with the same code as in Section 3.3, with k_g = 0.863 W/(mK), (ρc)_g = 3.8 MJ/(m³K), k_s = 3.18 W/(mK), and α_s = 10⁻⁶ m²/s. The values of θ ₁₀₀ and of θ _{100_approx} are reported on the right vertical axis. The figure shows that the approximation obtained by Equation (6) is very good, even for rather short times.

3.6. Estimate of α_s with Non-Uniform Temperature Rise of the BHE Surface

We repeated the estimate of α_s, with the same solver, absolute tolerance and mesh as in Section 3.2, by imposing, on the BHE surface, the time-dependent distribution of temperature rise obtained by applying Equation (6) and the experimental values of θ_s1(t). The lowest RMSD between the values of θ_s2(t), θ_s3(t), θ_s4(t) resulting from simulations and the experimental ones was obtained with α_s = 1.05 × 10⁻⁶ m²/s, and was equal to 0.0190 °C. The simulation with this value of α_s was repeated with a regularly refined mesh with 113,280 triangular elements, yielding values of θ_s2(t), θ_s3(t), and θ_s4(t) nearly identical to those obtained with the previous mesh. Due to the very high similarity with the conditions considered in Section 3.2, further checks of mesh independence did not seem necessary. Therefore, the result of the estimate of the thermal diffusivity of the sand is α_s = 1.05 × 10⁻⁶ m²/s, and the RMSD between the simulated values and the experimental ones is 0.019 °C. The agreement between the experimental and the simulated values of θ_s2(t), θ_s3(t), θ_s4(t) is illustrated in Figure 9. The figure also shows that the adiabatic condition at the boundary of the sandbox, employed in our simulations, reproduces the irregular experimental curve of θ_s5(t) better than the condition of constant temperature θ_s5 = 0. We also performed simulations with the boundary condition of uniform time-dependent temperature change of the lateral surface, equal to the experimental values of θ_s5(t), but we obtained a higher RMSD between the simulated and the experimental values of θ_s2(t), θ_s3(t) and θ_s4(t), with respect to the simulations with the adiabatic boundary. Therefore, we adopted the boundary condition of the adiabatic lateral surface and we recommend this choice for the validations of simulation codes by comparison with the results of the sandbox experiment.

3.7. Estimate of k_s with the Values of k_g and α_s Determined Previously

We repeated the estimate of k_s with the same method, solver, absolute tolerance, mesh and value of (ρc)_g as in Section 3.3, with the value of k_g determined in Section 3.4 and the value of α_s determined in Section 3.6, namely k_g = 0.863 W/(mK) and α_s = 1.05 × 10⁻⁶ m²/s. The lowest RMSD between the values of θ₁(t), θ₂(t), θ_s1(t), θ_s2(t), θ_s3(t), θ_s4(t) resulting from simulations and the experimental ones were obtained with k_s = 3.22 W/(mK), and was equal to 0.131 °C. The simulation with k_s = 3.22 W/(mK) was repeated with a mesh obtained by a regular refinement of the previous one, with 208,608 triangular elements. The refined mesh yielded nearly identical results. The RMSD between the results obtained by the different meshes was 0.00011 °C. Due to the very high similarity with the conditions considered in Section 3.3, further checks of mesh independence did not seem necessary. Therefore, the result of the estimate of the thermal conductivity of the sand is k_s = 3.22 W/(mK), which corresponds to (ρ c)_s = 3.0667 MJ/(m³K), and the RMSD between the simulated values of θ₁(t), θ₂(t), θ_s1(t), θ_s2(t), θ_s3(t), θ_s4(t) and the experimental ones is 0.131 °C. A comparison between the experimental and the simulated values of θ₁(t), θ₂(t), θ_s1(t), θ_s2(t), θ_s3(t), and θ_s4(t) is illustrated in Figure 10. The figure shows a very good agreement with the experimental values, except for a slight overestimation of θ₁(t), θ₂(t), and θ_s1(t) between 20 and 600 min, and a slight underestimation of θ_s1(t) for t > 1800 min.

3.8. Estimate of (ρ c)_g

The overestimation of θ₁ and θ₂ between 20 and 600 min illustrated in Figure 10 disappears in the long time. Therefore, it can be due to either an underestimation of the heat capacity of water, i.e., of the additive volume of water considered to take into account the thermal inertia of the external circuit, or to an underestimation of the heat capacity of the grout. Since the time evolutions of θ ₁ and of θ ₂ during the first 20 min are reproduced correctly, the only possibility remains an underestimation of the volumetric heat capacity of the grout. Values of this quantity higher than that recommended by Beier [30], namely 3.8 MJ/(m³K), are unusual. However, there are unavoidable differences between the experimental reality and the scheme thereof adopted in simulations, such as the difference between the real thermal insulation placed at the end surfaces of the sandbox, not described in Reference [35], and the condition of adiabatic surfaces. These differences could require values of some parameters slightly different from the real ones in order to reproduce the experimental results. Since there were no other ways to eliminate the overestimation of θ ₁ and of θ ₂ between 20 min and 10 h without reducing the overall accuracy, we determined the value of (ρ c)_g that minimizes the root mean square deviation between the simulation results and the experimental data without changing the parameters already optimized. We employed the same method, solver, absolute tolerance and mesh as in Section 3.7, and the following values of the thermal properties of grout and sand: k_g = 0.863 W/(mK), k_s = 3.22 W/(mK), and (ρ c)_s = 3.0667 MJ/(m³K). The result was (ρ c)_g = 4.6 MJ/(m³K). A color map of the temperature distribution at the final instant of time, with isothermal lines, is presented in Figure 11. The nearly perfect agreement between the experimental time evolutions of θ₁, θ ₂, θ _s1, θ _s2, θ _s3, θ _s4 and those obtained by the 2D simulation is shown in Figure 12. The root means square deviations between simulated and experimental values for θ ₁ and θ ₂ are 0.058 and 0.053 °C, respectively. The root mean square deviation for the complete set of values of θ ₁, θ ₂, θ _s1, θ _s2, θ _s3, θ _s4 is 0.082 °C. The accuracy obtained in reproducing the time evolution of θ ₁ and θ ₂ even during the first hour is evidenced in Figure 13.

4. Validation of the Simulation Codes

The excellent agreement between the experimental time evolutions of θ ₁, θ ₂, θ _s1, θ _s2, θ _s3, θ_s4 and the simulated ones is a clear signal of the accuracy of our 2D simulation codes. Nevertheless, we performed a validation of the simulation code employed in Section 3.3, Section 3.7 and Section 3.8 to reproduce the time evolutions of θ ₁, θ ₂, θ _s1, θ _s2, θ _s3, and θ_s4 by a comparison with the outcomes of the analytical BHE model by Man et al. [27]. In this model, the BHE is represented by a solid cylinder with the same thermal properties as the ground, which contains a heat-generating surface placed at a suitable radial distance r₀ from the BHE axis. The mean temperature of the surface at r = r₀ represents that of the BHE fluid. In order to improve the accuracy of the model during the first hours, where the real BHE geometry and heat capacity have a relevant effect on the mean-fluid-temperature rise, θ _fm, we considered for our validation a BHE with the geometry of the sandbox experiment and the thermal properties of the BHE materials equal to those of the ground, selected as k_s = 2.82 W/(mK), (ρ c)_s = 3.2 MJ/(m³K), and thus α_s = 8.8125 × 10⁻⁷ m²/s. The finite element simulations of the real BHE were performed by the same method, solver, absolute tolerance and mesh as in Section 3.3, Section 3.7 and Section 3.8.

The analytical expression of the fluid temperature rise, θ, as a function of the radial coordinate, r, and time, t, determined by Man et al. [27], is

$$\theta \left(r,t\right)=-\frac{q_b}{4{\pi }^2{k}_s}{\displaystyle {\int }_0^{\pi }Ei\left\{-\frac{r^2+r_0^2-2\hspace{0.17em}r\hspace{0.17em}{r}_0\cos \phi }{4\hspace{0.17em}{\alpha }_{s}t}\right\}}\hspace{0.17em}d\phi$$

(7)

where Ei is the exponential integral function and φ is an angular coordinate. The value of r₀ that minimizes the root mean square deviation between the values of θ _fm evaluated by the 2D finite-element simulation and those determined by the numerical integration of Equation (7) is 2.37 cm. With this value of r₀, we evaluated, by Equation (7), the time evolutions of θ _fm, of the mean temperature rise of the external surface of the BHE, θ _b, and of the mean temperature rise of the surface with r = 30.5 cm, θ _30.5, and we compared the results with those obtained by the 2D finite-element simulation. The comparison is illustrated in Figure 14, where the results of the 2D simulation are denoted by the subscript sim and those obtained by Equation (7) are denoted by the subscript Man. In the figure, the curves obtained by the finite-element simulation appear as coincidental with those obtained through the analytical solution. However, some discrepancy occurs during the first two hours, especially for θ _fm, due to the non-perfect accuracy of the analytical model for short times, as is evidenced in Figure 15, which refers to the time interval from 0 to 120 min.

The RMSD between the results of the finite-element simulation and those obtained by the analytical solution, in the whole time interval from 0 to 3106 min, is 0.016 °C for θ _fm, 0.005 °C for θ _b, and 0.004 °C for θ _30.5. In the time interval from 120 to 3106 min, the RMSD is 0.004 °C for θ _fm, 0.003 °C for θ _b, and 0.004 °C for θ _30.5.

Thus, the validation confirms the excellent accuracy of our 2D simulation code.

5. Validation of the Estimate by 3D Simulations

The results of our estimate, namely k_g = 0.863 W/(mK), (ρ c)_g = 4.6 MJ/(m³K), k_s = 3.22 W/(mK), (ρ c)_s = 3.0667 MJ/(m³K), and total volume of water in the system equal to that contained into the BHE tubes plus 8 L, were checked by 3D simulations of the system.

We considered, as a computational domain, a rectangular parallelepiped with square bases, having a length of 18.32 m and bases with sides of 1.83 m. All the external surfaces of the domain were assumed as adiabatic, except for the inlet and the outlet cross section of each pipe. The water velocity was set as uniform, with a value of 0.3358 m/s, and only the local energy balance equation was solved in the fluid domain. The time dependent inlet temperature was imposed through the energy balance on a thermally insulated tank with a volume of 8 L that receives a constant thermal power $\dot{Q}$ = 1051.6 W. Under the assumption that the rate of change of T_out can be neglected as long as the difference T_in–T_out is significantly lower than its asymptotic value, the energy balance equation yields [53]

$$T_{in}\left(\tau \right)=T_{out}\left(\tau \right)+\frac{\dot{Q}}{\dot{m}{c}_{pw}}\left(1-e^{-\frac{\dot{m}{c}_{pw}}{C_w}\tau }\right)$$

(8)

where T_out is the surface average of the fluid temperature at the outlet section, $\dot{m}$ is the mass flow rate, 0.19614 kg/s, c_pw is the specific heat capacity at a constant pressure of water, 4.1798 kJ/(kg K)), and C_w is the heat capacity of the water contained in the tank, equal to 33.293 kJ/K. A heat flux per unit area given by the product of the convection coefficient and the temperature difference between fluid and solid surface was imposed at the fluid-solid interface.

The initial temperature in the solid domain was set uniform and equal to 22.068 °C, i.e., to the average between the mean fluid temperature and the temperatures measured by the sensors S1, S2, S3, S4, and S5. The initial temperature in the fluid domain was set as a linear function of the distance from the inlet section, with extreme values equal to T_in and T_out, i.e., 22.211 and 21.978 °C.

In analogy to other papers where 3D simulations of BHEs with real geometry were performed [53,54,55,56], the coordinate z along the BHE length was replaced by a rescaled one, z_r, in order to obtain a more compact shape of the computational domain. After some preliminary simulations, we selected the rescaling factor 40, i.e., the relation z_r = z/40. As a consequence, we defined a reduced velocity equal to the real velocity divided by 40, and an anisotropic thermal conductivity, for all the media, with the axial component equal to the other components divided by 1600 [53,54,55,56].

Finite-element 3D simulations with k_g = 0.863 W/(mK), (ρ c)_g = 4.6 MJ/(m³K), k_s = 3.22 W/(mK), (ρ c)_s = 3.0667 MJ/(m³K), and 8 L of external water were performed through COMSOL Multiphysics, with the solver PARDISO, absolute tolerance 10⁻⁴, and two unstructured meshes: Mesh 1, having 661,872 tetrahedral elements, and Mesh 2, having 1,683,706 tetrahedral elements. An illustration of the system, with Mesh 2, is reported in Figure 16a, and a particular of the upper surface, with different colors for different materials, is represented in Figure 16b.

The time evolutions of T_in, T_out, T_s1 and T_s2 obtained by Mesh 1 and by Mesh 2 are compared withthe experimental time evolutions in Figure 17, where Mesh 1 and Mesh 2 are denoted by M1 and M2. The figure shows that increasing the number of elements yields a reduction in the discrepancy between the simulation results and the experimental ones, and that the time evolutions obtained by Mesh 2 are in fair agreement with the experimental ones. The root means square deviations between the simulated and the experimental values of T_in and of T_out are 0.07 °C and 0.06 °C for Mesh 2, both 0.14 °C for Mesh 1. The overall root mean square deviation from the experimental values is 0.08 °C for Mesh 2 and 0.12 °C for Mesh 1. Therefore, the results of the 3D simulations confirm the accuracy of our estimate.

6. Conclusions

By means of 2D finite-element simulations, we have provided a new estimate of the thermal properties of the grout and of the sand of the laboratory BHE model studied experimentally by Beier et al. [35]. In our analysis, we have considered all the available details of the experiment, including the presence of a thin aluminum pipe at the BHE boundary and the difference between the local temperatures of the sand in correspondence to the axis of the inlet pipe, measured in Reference [35], and the mean temperatures at the same radial distances. All the experimental data have been kindly provided to us by the first author of Reference [35].

First, we performed a preliminary estimate of the thermal conductivity of the grout, k_g, by a stationary simulation, with thermal conductivity of the sand k_s = 2.82 W/(mK), as stated in Reference [35]. Then, we performed a preliminary estimate of the thermal diffusivity of the sand, α_s, by a simulation of the transient heat transfer in the sand, by assuming a uniform temperature of the external BHE surface. With the obtained values of k_g and α_s, we performed a preliminary estimate of the thermal conductivity of the sand, k_s, under the assumption that the volumetric heat capacity of the grout is (ρ c)_g = 3.8 MJ/(m³K), as estimated by Beier [30]. With the estimated value of k_s, we repeated the estimate of k_g. Then, we repeated the estimate of α_s by considering the non-uniform time-dependent temperature distribution at the BHE surface, and we repeated the estimate of k_s with the new value of α_s. Finally, we estimated the value of (ρ c)_g that minimizes the root mean square deviation between the measured temperatures and the simulated temperatures, without changing the other parameters. The principal 2D simulation code has been validated by comparison to the analytical model of Man et al. [27]. The accuracy of the estimate has been confirmed by 3D simulations.

The results of our estimate are: heat flux per unit length supplied to the BHE equal to 57.401 W/m, k_g = 0.863 W/(mK), (ρ c)_g = 4.6 MJ/(m³K), k_s = 3.22 W/(mK), α_s = 1.05 × 10⁻⁶ m²/s, heat capacity of the circuit external to the BHE equivalent to that of eight additional liters of water in the BHE. These parameters, applied to the geometry of the BHE and of the temperature sensors defined in Reference [35], have allowed us to reproduce the time evolution of the mean temperature of water in the inlet tube, of water in the outlet tube, of the sensors at the BHE external surface, and of the sensors at the radial distances 30.5 and 51.5 cm from the BHE axis with an overall root mean square deviation of 0.082 °C. The root mean square deviations obtained for the mean water temperature in the inlet pipe and that in the outlet pipe are 0.058 °C and 0.053 °C, respectively. The set of values reported above is therefore recommended for future validations of simulation codes by comparison withthe sandbox experiment [35].

The higher thermal conductivity of the grout, with respect to that measured in Reference [35], could be due to absorption of water from the wet sand. The very high value of the volumetric heat capacity of the grout found in this paper could be due, in part, to heat losses from the thermally insulated end surfaces of the sandbox, considered as adiabatic in simulations.