An improved method for upscaling borehole thermal energy storage using inverse finite element modelling

Accepted Manuscript An improved method for upscaling borehole thermal energy storage using inverse finite element modelling K.W. Tordrup, S.E. Poulsen, H. Bjørn PII:

S0960-1481(16)31060-6

DOI:

10.1016/j.renene.2016.12.011

Reference:

RENE 8348

To appear in:

Renewable Energy

Received Date: 21 December 2015 Revised Date:

2 December 2016

Accepted Date: 5 December 2016

Please cite this article as: Tordrup KW, Poulsen SE, Bjørn H, An improved method for upscaling borehole thermal energy storage using inverse finite element modelling, Renewable Energy (2017), doi: 10.1016/j.renene.2016.12.011. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

3

K. W. Tordrup*, S. E. Poulsen and H. Bjørn.

4

Centre of Applied Research and Development in Building, Energy & Environment, VIA University College,

5

8700 Horsens, Denmark

6

*

7 8

Abstract

9

natural variability of thermal conductivity and heat capacity in the storage volume. We present an improved

10

method for upscaling a pilot BTES to full scale and apply the method to an operational storage in Brædstrup,

11

Denmark. The procedure utilizes inverse 3D finite element method (FEM) modelling of distributed

12

temperature measurements inside the BTES for inferring the thermal properties of the subsurface. We find

13

that individual geological layers can be distinguished in terms of their heat capacities and thermal

14

conductivities using inverse modelling. The depth integrated estimate of thermal conductivity differs

15

significantly from that obtained from a single thermal response test (TRT) at the site. As such, we find

16

significant scaling effects in terms of the subsurface thermal conductivity distribution which are expected to

17

be further amplified in an expansion of the pilot BTES to full scale. The methodology presented in this paper

18

therefore provides an improved basis for upscaling pilot BTES systems. The operational data and BTES

19

temperature measurements are published with the present paper in the supplementary material.

20

Keywords: borehole thermal energy storage, numerical modelling, inverse modelling, model validation,

21

thermal conductivity, dimensioning

Corresponding author: [email protected] (K. W. Tordrup)

RI PT

2

An improved method for upscaling borehole thermal energy storage using inverse finite element modelling

TE D

M AN U

SC

Dimensioning of large-scale borehole thermal energy storage (BTES) is inherently uncertain due to the

EP

1

1. Introduction

Large-scale thermal solar collector arrays are currently being integrated into the district heating networks in

24

Denmark in order to increase the fraction of renewable energy of total heat consumption. Since solar energy

25

production peaks in the summer when heat consumption is low, storage is required for effective balancing of

26

seasonal fluctuations in energy demand and consumption. In addition, the increased production of wind

27

power makes it economically feasible to produce heat by means of heat pumps or electrical boilers when the

28

price of electricity is low, which further increases the need for heat load balancing. In Denmark, district

29

heating networks with large solar collector arrays are typically balanced by a combination of accumulation

30

tanks and hot water pit storage. Currently, the district heating networks in Denmark incorporate more than

31

80, large-scale solar thermal arrays with either pit or tank storage systems [1].

AC C

22 23

1

ACCEPTED MANUSCRIPT Borehole thermal energy storage is a viable alternative to hot water pit storage. For a review of BTES

33

technology we refer to [2] and [3]. Practical experiences from existing BTES systems demonstrate different

34

levels of success as outlined in the following. During the second year of operation of the Anneberg BTES in

35

Sweden [4], electricity-based heat consumption was reduced by just one-third of the long-term target of 70%.

36

Moreover, practical problems related to the construction of the BTES were encountered including leaky heat

37

exchanger pipes and removal of contaminated soil from the storage after the system had been taken into use.

38

In Attenkirchen, Germany the BTES technology has been combined with pit storage in order to improve the

39

responsiveness of the storage [5]. In Crailsheim, Germany a conventional BTES has demonstrated a solar

40

fraction of 35.8% after four years of operation, which is well below the design target of 51% [6,7]. In Torino

41

Italy, an experimental BTES with 4 borehole heat exchangers (BHEs) was constructed and put into operation

42

in 2014 [2]. The authors find that during the first year of operation the solar collectors were able to produce

43

just two thirds of the forecasted heat production. Correspondingly, during discharge of the BTES 17-22% of

44

the injected energy was reclaimed. The Drakes Landing BTES in Canada stores solar energy in

45

approximately 34,000 m3 of soil and supplies 52 energy efficient homes with district heating [8]. The Drakes

46

Landing BTES facility has demonstrated solar fractions of 97% after 5 years of operation, which is well

47

above the target of 90%.

48

The local district heat and power company in Brædstrup, Denmark has constructed a 19,000 m3 pilot BTES

49

in order to store excess heat production from an 18,600 m2 thermal solar collector array [9]. The aim of the

50

pilot project is to evaluate the ability of the BTES for load balancing of the district heating network. The

51

BTES is scheduled to expand the storage volume to approximately 200,000 m3, given that it performs

52

satisfactorily. The initial, thermal characterization and predicted capacity of the Brædstrup BTES are based

53

on a single TRT-based soil thermal conductivity estimate. Since the soil thermal conductivity quantifies the

54

ease with which heat can be abstracted from and dissipated in the BTES it represents a vital performance

55

parameter. The volume of soil sampled in a TRT is exceedingly small relative to the total BTES volume and

56

the associated soil thermal conductivity estimate may prove to be unrepresentative for the total storage

57

volume. Such scaling effects potentially bias model forecasts of BTES thermal behavior and storage

58

capacity. Accurate prediction of heat injection and abstraction rates is vital for ground source heating

59

systems due to the relatively high costs of the technology [10]. The subsurface characterization required to

60

make such predictions must encompass the thermal properties of the storage medium and boreholes as well

61

as the hydrogeological conditions in the presence of groundwater flow [11,12].

62

In this study, we present operational data and storage temperature measurements from the pilot BTES in

63

Brædstrup, Denmark in an improved inverse modelling based characterization of the thermal properties of

64

the BTES. The characterization serves to improve the basis for a potential future expansion of the BTES to

AC C

EP

TE D

M AN U

SC

RI PT

32

2

ACCEPTED MANUSCRIPT full scale. Inverse modelling of measured BTES temperatures yielding improved thermal conductivity

66

estimates has not been demonstrated in previous literature.

67

The paper is organized as follows. In Section 2 we present the Brædstrup BTES and its integration into the

68

local combined heat and power plant. In Section 3 we describe the preprocessing of operational data from the

69

first 510 days of operation. In Section 4 we develop a 3D numerical finite element model for transient

70

simulation of BTES soil temperatures. The BTES model is calibrated using temperature observations in five

71

boreholes, at twenty different levels, placed inside and in the vicinity of the storage. The thermal

72

conductivity and heat capacity of six identified lithological units are estimated in the calibration. In Section 5

73

we present calibration results and validate the model by comparing predicted and observed BTES

74

temperatures outside the calibration dataset. To investigate potential scaling effects, the calibrated model

75

estimate of the depth integrated soil thermal conductivity is compared to the corresponding estimate obtained

76

from a single thermal response test. The predicted BTES energy balance calculated with the calibrated six-

77

layer model is compared to the corresponding prediction assuming a homogeneous subsurface thermal

78

conductivity equal to that obtained from the thermal response test. We explore the degree to which variations

79

in the subsurface thermal conductivity impact scaling of the BTES to full size, and how this affects model

80

projections of the ability of the BTES to store and release heat when load balancing is required. Conclusions

81

are drawn in Section 6.

2. System description

M AN U

SC

RI PT

65

Brædstrup District Heating delivers heat to approximately 1,500 households. Annually this amounts to 45

84

GWh of heat. The heat is distributed to consumers at a temperature of 72-80°C depending on the season. The

85

plant has two gas boilers as well as two gas generators that both run on natural gas. As a pilot project a

86

19,000 m3 BTES has been constructed for seasonal storage of surplus heat from an 18,600 m2 solar collector

87

field. Short term load balancing is supplied by two steel accumulation tanks with a combined volume of

88

7,500 m3 and heat is extracted from the storage using a 1.2 MW heat pump (see Fig. 1).

AC C

EP

TE D

82 83

3

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

89

Figure 1: Simplified diagram of the Brædstrup combined heat and power plant.

91

In addition, the facility has a 10 MW electric boiler, which can consume surplus power from the electrical

92

grid. The pilot project facilitates approximately 20 % solar coverage of the total heat demand on an annual

93

basis. During days with high solar production, the solar collectors first cover the consumer demand.

94

Secondly, the water-based accumulation tanks are charged and any excess production is subsequently stored

95

in the BTES. The BTES is discharged by the heat pump, that boosts the temperature to 85 °C in one

96

compression cycle. The lower cut-off temperature from the BTES is 16 °C in order to maintain an acceptable

97

coefficient of performance.

98

Depending on the outcome of the pilot project, a full scale plant with 50,000 m2 thermal solar panels and

99

200,000 m3 BTES is to be constructed. The ultimate goal is 50 % solar coverage of the total annual heat

EP

TE D

90

demand.

101

2.1 Geological and geothermal setting

102

A soil profile for the BTES has been established from two drillings (see Figure 2).

4

AC C

100

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

103 104

Figure 2: Lithological profile and placement of the BHEs (grey) and a temperature monitoring drilling (red) in the BTES.

105

The BHEs are covered by a 0.5 m insulating layer of mussel shells and 0.5 m of soil fill. The BHEs extend to

106

45 m beneath the insulation and separate boreholes for measuring soil temperature extend to 59 m beneath

107

the insulation. The drillings show sediments that consist of glacial deposits below a thin layer of topsoil. Soft

108

strata of clay till and silt are encountered down to approximately 2-3 m depth. The clay till and silt are

109

underlain by a strata of clay till down to 8–10 m depth at which deposits of sand and gravel are found. At

5

ACCEPTED MANUSCRIPT 19–20 m depth a layer of clay till extends down to 23.5–28 m depth. The clay till is underlain by deposits

111

alternating between fine sand and sandy silt with embedded layers of clay till and sand till. Below 41 m

112

depth, alternating layers of sand and silt are found. The sediments are primarily dry between 3-50 m depth

113

with 2-8% water content. Temporary seepage of groundwater was observed in the upper clay till during the

114

drilling due to infiltration of melting snow. During the drilling, groundwater was encountered at 49.3 m

115

depth in a standpipe that reaches a depth of 51 m. However, at a later sounding no groundwater table was

116

encountered. The regional groundwater table is observed at 69 m depth in a 102 m deep well situated 20 m

117

from the BTES [13]. We conclude that the alternating sand and silt comprise a perched aquifer in which a

118

temporary water table forms when infiltration is high.

119

Two weeks following drilling completion, a thermal response test of a single BHE was carried out. The

120

standard line-source based interpretation of the test implies an effective soil thermal conductivity along the

121

length of the BHE of 1.42 W/m/K and a borehole thermal resistance of 0.172 K·m/W. The volumetric heat

122

capacity was estimated to 1.9 MJ/m3/K from standard table values of unsorted sand with water contents

123

similar to that found in the BTES volume.

124

2.2 BTES description

125

The BTES consists of 48 boreholes arranged in a hexagonal pattern with a distance of 3 m between the

126

boreholes (see Figure 3), each of which contains a 2U heat exchanger.

AC C

EP

TE D

M AN U

SC

RI PT

110

127 128

Figure 3: The BTES layout including five temperature probe boreholes T1-T5 and the central manifold (black disc).

6

ACCEPTED MANUSCRIPT For all temperature probes T1-T5 the temperature is recorded at z = 0, 1, 2, 3, 4, 9, 14, 19, 24, 29, 34, 39, 44,

130

48, 49, 50, 51, 52, 53 and 59 m below the bottom of the insulation.

131

The 96 U-tube heat exchangers in the 48 boreholes are arranged in 16 parallel strings with 6 heat exchangers

132

in each, as indicated in Figure 4.

AC C

EP

TE D

M AN U

SC

RI PT

129

133 134

Figure 4: Configuration of the 16 BHE array strings.

7

ACCEPTED MANUSCRIPT Thus, the two U-tube heat exchangers in each borehole are connected to different strings. This way, if a

136

string is taken out of operation, the boreholes in that string will still be in use. The BTES is charged from the

137

center of the storage, and during discharge the direction of the flow is reversed such that discharging

138

proceeds from the outside of the storage volume. Since the tubes are located beneath the insulation there is

139

no risk of freezing and therefore the heat carrier fluid is pure water to avoid soil contamination in the event

140

of a leak.

141

Soil temperatures are recorded using PT100 sensors at five locations as indicated in Figure 3. Additionally,

142

the inlet and outlet temperatures of the 16 BHE arrays, as well as the total BTES flowrate are recorded. The

143

flow is equally distributed between the 16 strings. All measurements have been recorded for 510 days

144

starting on May 22nd 2012.

SC

RI PT

135

3. Data preprocessing

The temperature time series at all positions shown in Figure 3 as well as inlet and outlet temperatures for the

147

16 BHE array strings and the total flowrate to the BTES have been extracted from the SCADA system at the

148

district heating facility in Brædstrup. The measurements were recorded as hourly average values, however

149

due to an error in the SCADA system, the data for the first 11 days were stored as 5 minute average values.

150

In order to eliminate spurious fluctuations and minimize the computational time for the FEM model, all

151

temperature time series are aggregated to daily averages. Since the thermal response of the soil occurs on a

152

timescale of days the effect of this is minimal. In order to ensure conservation of energy, the flowrate to the

153

BTES is aggregated to daily values as indicated in Eq. (1)

154


 = ∑

155

where
is the volumetric flowrate [m3/h]. Subscripts h and l refer to the raw high resolution data and low

156

resolution aggregated values respectively. ∆, and ∆, are the difference between the inlet and outlet

157

temperature [K] for the i’th BHE array string in high and low resolution, respectively, and ∆ = 24 ℎ

158

specifies the resolution of the aggregated data.

159

An automatic procedure has been implemented to remove outliers from the dataset. For each of the

160

aggregated time series the first derivative is evaluated numerically. A data point is deemed an outlier if the

161

derivative departs from its average value by more than two standard deviations. In these cases the outlier is

162

replaced by an interpolated value.

163

Finally, all the time series have been inspected manually. For the measurements of the BHE inlet and outlet

164

temperatures it was observed that the outermost temperature sensor in string 11 (see Figure 4) recorded a

165

constant value of 150°C for a period of 450 days. This measurement has been replaced by the average of the

TE D

M AN U

145 146

(1)

AC C

EP

∆   ∑
 ∆, ,  ∆ ∆

,



8

ACCEPTED MANUSCRIPT outermost temperatures of the remaining 15 strings. Since the inlet and outlet temperatures are all within a

167

few degrees of each other this will not significantly impact the model predictions.

168

The aggregated temperature profiles for the five temperature probe boreholes T1-T5 at selected depths and

169

aggregated total flow are plotted in Figure 5.

AC C

EP

TE D

M AN U

SC

RI PT

166

170

9

ACCEPTED MANUSCRIPT Figure 5: (a-e): Aggregated temperature profiles for T1-T5 at selected depths. (f): Aggregated total flow and inlet temperature to the storage.

173

During the summer, excess heat production from the solar collectors is stored in the BTES causing soil

174

temperatures to increase to in excess of 50°C. During winter, the BTES is discharged and temperatures drop

175

to slightly above initial levels. Below the storage, the thermal disturbances are moderate and phase-shifted

176

with respect to that of higher layers. T4 is placed approximately 11 m from the outermost BHEs and as

177

expected the temperature disturbance is dominated by the seasonal cycle of the solar insolation in the upper

178

layers. Effects from storage in the BTES are evident in lower layers as a slight increase in soil temperature

179

(see Fig. 5d). The unprocessed operational data and measured storage temperatures are available in the

180

supplementary material.

181

Of the 100 soil temperature time series, 14 have been omitted from the study. All temperature measurements

182

at z = 0 m are omitted since they reside on the model boundary. It was discovered during the analysis that the

183

recorded temperatures at T2 at z = 48 m, are exactly identical to those recorded at T2 at z = 52 m. It was

184

judged that the measurements were most likely to be from z = 52 m, so the measurements at z = 48 m were

185

discarded. In addition, it appeared that the measurements at T2 at z = 34 m and T5 at z = 44 m have been

186

interchanged in the SCADA system. Since it was not possible to achieve definite confirmation of this

187

suspicion, both measurements were discarded as a precautionary measure. The time series at T4 at z = 34 m

188

registered a constant value of -20°C for 243 days and was therefore discarded. Since T4 is located outside the

189

storage, the uppermost temperature sensors are sensitive to solar insolation which is not simulated by the

190

FEM model. As such, the measurements at z = 1, 2, 3 and 4 m at T4 are not included in the calibration.

191

During the initial calibration attempt it was found that the time series at T5 at z = 1 m caused problems with

192

convergence. Presumably since this measurement is located at the very center of the storage and only 1 m

193

below the insulation our model does not adequately capture the thermal dynamics at this point, and so it is

194

discarded from the calibration. Finally, at T3 the temperatures recorded at levels z = 24 m, 29 m, 34 m, and

195

39 m stick at a plateau for 15 days during the second heating season at t = 430–444 d. As these data points

196

are used solely to validate the model and are not included in the calibration, they are simply replaced by

197

interpolated values.

AC C

EP

TE D

M AN U

SC

RI PT

171 172

4. Methods

198 199

4.1 Finite element model

200

The finite element software FEFLOW [14] is utilized for calculating the subsurface temperature response in

201

and near the BTES. FEFLOW uses the finite element method to solve the governing equation for transient

202

thermal conduction in three dimensions:

203

ρc

10

"# ! "$

= ∇ ∙ '∇T + Q

(2)

ACCEPTED MANUSCRIPT where (ρc)b is the bulk volumetric heat capacity [J/m3/K], T the temperature [K], t the time [s], and λ is the

205

bulk thermal conductivity tensor [W/m/K]. Q is the heat generation rate [W/m3], which in the present case is

206

due to the BHEs. The BHEs are modelled using the quasi-stationary analytical model developed by Eskilson

207

and Claesson [15] using the physical data indicated in Figure 6 and a grout thermal conductivity of λg = 1.44

208

W/m/K. The quasi-stationary BHE model is restricted by the line source time criterion as well as the time

209

required to circulate fluid through the pipes [15]. The full FEM model uses 536,130 elements with 7659

210

elements in each of the 70 model layers, and the vertical discretization is 1 m.

TE D

M AN U

SC

RI PT

204

211

Figure 6: Cross section of a borehole heat exchanger. The outer diameter of the U-tubes is indicated and the wall thickness is 2.9 mm.

214

The horizontal distance from the BTES to the model boundary is determined from a line source based

215

estimate of the radial temperature disturbance during the simulated period, such that the modelled

216

temperature change at the boundaries of the model is negligible. The temperature change at the boundary of

217

the model is below 0.3 °C in all simulations carried out in this study. The model extends 20 m horizontally

218

from the outermost BHEs and vertically from immediately below the insulation to 25 m below the BTES.

219

The node density in the central region is roughly 3.8 times that of the surrounding area for increased

220

accuracy in the storage region (see Fig. 7). The thicknesses of the geological layers in the model are based on

221

the lithological profile in Figure 2. We use a triangular prism mesh, and the six nearest horizontal nodes to a

222

BHE-node are chosen at a distance of 6.13·D/2 for optimum accuracy [6]. The second-order Adam-

223

Bashforth predictor-corrector scheme is utilized for automatic time-stepping.

AC C

EP

212 213

11

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

224

Figure 7: The FEM model mesh. The extent of the geological layers is indicated by alternating shades of grey. BHE positions are indicated in red.

227

The error tolerance for the BHE model is set to 10-5 and the maximum allowed number of iterations is set

228

equal to 1. Finally, since FEFLOW cannot handle reversing the flow direction through a string of BHEs

229

under discharging, a plugin has been developed with the FEFLOW API (Application Programming

230

Interface) to explicitly link the boreholes in the correct order under charging and discharging, respectively.

231

4.1.1 Initial- and boundary conditions

232

The undisturbed ground temperature was measured to 7.9-8.2 °C in the depth interval penetrated by the

233

temperature probes, hence the initial temperature is set to 8 °C throughout the model volume. The thermal

234

gradient is set to zero along all model boundaries

235

∇T = 0,

236

implying that heat transport across the boundary is neglected. In particular, this implies that solar insolation

237

and heat flow from the Earth’s interior are not simulated in the model. Since the storage is insulated from

238

above, the effect of solar irradiance is negligible compared to the direct heat injection into and abstraction

239

from the BTES. In addition, any upward heat loss through the insulation is not included in the model. While

240

the effect on the overall energy balance will be negligible if the insulation is good, we cannot expect the

241

model to provide good predictions of the temperatures at the upper model boundary (z = 0). The heat flow

242

from the Earth’s interior is approximately 35 mW/m2 [16] yielding approximately 0.1 MWh/year for the

243

BTES storage volume, which is insignificant relative to the stored and abstracted energy during operation.

12

AC C

EP

TE D

225 226

(3)

ACCEPTED MANUSCRIPT 4.2 Parameter estimation

245

The bulk thermal conductivity and volumetric heat capacity are estimated for each of the six model layers

246

(see Fig. 7). Initial guesses for all layers are λ = 1.42 W/m/K as suggested by the thermal response test and

247

ρc = 1.9 MJ/m3/K as estimated from table values.

248

The parameter estimation is performed using the software PEST (Model-Independent Parameter Estimation

249

and Uncertainty Analysis) [17]. PEST utilizes the Gauss-Marquardt-Levenberg algorithm for minimizing the

250

objective function

251

56789 4 ∑. , = ∑ -,. ,. − ,.

252

where  / and  1 indicate measured and simulated temperatures respectively. :;<=>? is the number of

253

temperature probes, and : is the number of temporal observations in the dataset for each temperature probe.

254

-,. is the weight assigned to individual observations.

255

The first 310 days of subsurface temperature measurements serve as calibration data for estimating the

256

thermal conductivity and volumetric heat capacity of the six soil layers. Temperature observations recorded

257

for the first 20 days are assigned a weight of zero (-,. in Eq. (4)) in the calibration to allow initial transients

258

to decay without affecting the objective function. For the remaining 290 days, all temperature measurements

259

are assigned unit weights in the calibration. Weights should be chosen so that they are inversely proportional

260

to the standard deviation of the corresponding observation [17]. Since all temperatures are measured with

261

identical sensors and under identical conditions, it is reasonable to assign identical non-zero weights to all

262

observations included in the calibration.

263

Since FEFLOW does not permit separate configurations of the two U-tube heat exchangers in each borehole,

264

the calibration has been performed separately for the two configurations in Figure 4. Parameter estimates are

265

calculated as the arithmetic average of the corresponding individual estimates for the two configurations.

266 267

5. Results and discussion

268

Table 1 together with the final estimate of each model parameter. For the average parameters the worst case

269

error estimate from the two model configurations is reported.

/

1

2

,

(4)

SC

3

AC C

EP

TE D

M AN U

3

RI PT

244

The estimated soil parameters for the six model layers and for both BHE array configurations are given in

270 271

13

ACCEPTED MANUSCRIPT Table 1: Calibration estimates and linear 95% confidence intervals for the soil thermal conductivity λave and volumetric heat capacity (ρc)ave for each layer in the model and for each BHE configuration Figure 4.

Configuration 1 (ρc)1

λ1

(ρc)2

λ2

(ρc)ave

λave

[MJ/m3/K]

[W/m/K]

[MJ/m3/K]

[W/m/K]

[MJ/m3/K]

[W/m/K]

0 – 3m

2.04±0.030 2.26±0.023

2.01±0.033

2.27±0.048

3m – 9m

2.27±0.025 1.41±0.026

2.26±0.018

1.37±0.017

9m – 20m

1.68±0.007 1.59±0.015

1.81±0.005

1.67±0.016

20m

–

2.18±0.032 1.48±0.025

2.13±0.034

1.48±0.034

–

1.95±0.021 1.76±0.032

1.93±0.018

–

1.83±0.025 2.41±0.041

1.89±0.020

41m 41m 70m 274

2.03±0.033

2.27±0.048

2.27±0.025

1.39±0.026

1.75±0.007

1.63±0.016

2.15±0.034

1.48±0.034

1.75±0.022

1.94±0.021

1.75±0.032

2.47±0.024

1.86±0.025

2.44±0.041

M AN U

26m 26m

Average

RI PT

Layer

Configuration 2

SC

272 273

For the upper clay till/silt at 0-3 m depth, the thermal conductivity estimate is relatively high compared to

276

values reported in the literature [18]. This is most likely due to the upper, perfectly insulating model

277

boundary condition. In the calibration process, the actual, upward heat loss is compensated by increasing the

278

thermal conductivity of the uppermost clay and till layer. Moreover, structural errors pertaining to

279

differences between the modelled and actual thickness of the thin uppermost clay till potentially bias the

280

thermal conductivity estimate. Thermal conductivity estimates for the geological layers in the depth interval

281

3-41 m depth correspond well to expectations [18]. The clay tills in the depth intervals 3-9 m and 20-26 m,

282

respectively, have lower thermal conductivities relative to the sand and gravel at 9-20 m and the fine sand at

283

26-41 m as expected. Conversely, the calibrated heat capacities are low for the sandy sediments and high for

284

the clay tills. The thermal conductivity estimate for the alternating sand and silt between 41 and 70 m depth

285

is relatively high. A drilling report from a 102 m deep well situated 20 m from the BTES, shows sand and

286

gravel at similar depth with an underlying layer of clay till. This observation indicates a certain degree of

287

inhomogeneity in the sediment composition and also explains the presence of the secondary water table

288

discussed in Section 2.1. The occurrence of periodic groundwater flow and the potential presence of coarser

289

grained materials may serve to explain the relatively high estimate of thermal conductivity.

AC C

EP

TE D

275

14

ACCEPTED MANUSCRIPT 5.1 Validation

291

In order to evaluate the predictive capabilities of the model, simulated soil temperatures were compared to

292

corresponding observations for the 200 days following the calibration period. Measurements at T1, T2, T3

293

and T5 between the surface and z = 44 m are included inside the storage volume. The distribution of the error

294

 / −  1 is indicated separately inside and outside the storage volume in Figure 8.

AC C

EP

TE D

M AN U

SC

RI PT

290

295 296

Figure 8: Histogram of prediction errors inside and outside the storage volume.

297

Inside the storage volume, the mean deviation of the errors is -0.76 °C and the variance is 1.5 (°C)2. The

298

negative bias in the simulation error indicates that soil temperatures are slightly overestimated inside the

299

storage volume, which is consistent with neglecting heat loss through the top insulation and thus trapping

300

heat in the model. Outside the storage volume the mean deviation is 0.66 °C and the variance is 0.16 (°C)2.

301

On average, soil temperatures outside the storage are slightly underestimated which may be an effect of

15

ACCEPTED MANUSCRIPT neglecting solar insolation in the area surrounding the BTES. In both cases, we note the bias is quite low and

303

as expected the variance is larger inside the storage volume than outside.

304

5.2 Scaling effects

305

The calibration of the 3D model yields improved estimates of soil thermal conductivity, relative to that of

306

thermal response testing, since a significantly greater soil volume is sampled in the full 3D FEM model. The

307

depth–integrated average of the thermal conductivity of the BTES volume is estimated as the thickness

308

weighted contribution of the individual layers, since heat transport mainly occurs parallel to the geological

309

layers. From the results in Table 1, the depth integrated estimate of the thermal conductivity of the storage

310

volume is 1.72 W/m/K which differs significantly from the value of 1.42 W/m/K obtained in the thermal

311

response test. In order to quantify the consequences of this discrepancy in terms of dimensioning the BTES

312

we compare the predictions of our calibrated model to a naive finite element simulation assuming a uniform

313

heat capacity of 1.9 MJ/m3/K and thermal conductivity of 1.42 W/m/K throughout the entire model volume.

314

The observed accumulated energy injection is compared to both models in Figure 9.

EP

TE D

M AN U

SC

RI PT

302

315

Figure 9: Accumulated energy injection calculated from observed fluid temperatures, calibrated model and naive model.

317

The calibrated model predicts the observed energy balance extremely well during the first heating season.

318

During the subsequent discharge period the model prediction slowly diverges from the observed accumulated

319

energy. During the first charging cycle heat is injected in a few pulses at a steady flowrate (see Fig. 5(f)).

320

Conversely, the initial discharge cycle proceeds with variable flow and in shorter pulses. Under these

321

conditions the assumptions of the quasi-stationary BHE model are violated and simulated fluid outlet

322

temperatures are accordingly inaccurate. The second heating cycle is again dominated by few injection

323

pulses with steady flowrates and the final error of the calibrated model is primarily an offset deriving from

324

the error accumulated during the initial discharge cycle. Conversely, the naive model fails to predict

325

accurately the observed accumulated energy during both charging and discharging. This is not due to the

AC C

316

16

ACCEPTED MANUSCRIPT limitations of the BHE model, but rather an indication that a homogenous model cannot encompass the full

327

thermal dynamics of the BTES. At the end of the second period of charging the observed injected energy is

328

equal to 616 MWh. The calibrated model predicts 591 MWh while the naive model predicts 539 MWh.

329

Thus, by including spatio-temporal temperature information in the calibration we are able to reduce the error

330

of the predicted energy balance from 12.5 % to 4.0 %.

331

Finally, we consider the predicted power injected into the BTES. In Figure 10 we compare the duration curve

332

of the daily average power predicted by the naive and calibrated models to the observed power.

M AN U

SC

RI PT

326

TE D

333

Figure 10: Duration curve for energy injection calculated from observed fluid temperatures and flow rates compared to model predictions.

336

As expected the naive TRT-based model consistently underestimates the rate at which energy can be injected

337

into the BTES since the thermal response test in the present case underestimates the thermal conductivity of

338

the subsurface as discussed previously.

339 340

6. Conclusions

341

procedure utilizes inverse 3D finite element modelling of distributed temperature measurements inside the

342

BTES for inferring the thermal properties of the subsurface. The forward model accepts measured flow rates

343

and fluid temperatures as input and computes outlet fluid and soil temperatures. We demonstrate that bulk

344

thermal conductivity and volumetric heat capacity of six identified stratigraphic units are well constrained by

345

non-linear regression utilizing subsurface BTES temperature measurements as calibration data. The error of

346

the energy balance prediction relative to observation is reduced from 12.5 % to 4.0 % when utilizing the

347

calibrated model relative to a naive TRT-based model parametrization. Consequently, the calibration of the

348

3D model yields improved estimates of the subsurface thermal conductivity distribution as well as storage

349

capacity compared to estimates based on a single TRT. In the present case, the depth-integrated thermal

EP

334 335

AC C

We develop an improved method for upscaling a pilot BTES in Brædstrup, Denmark to full scale. The

17

ACCEPTED MANUSCRIPT conductivity estimate calibrated with the full 3D model is 1.72 W/m/K compared to the TRT-based estimate

351

of 1.42 W/m/K. This implies improved predictions of the power with which the BTES may be charged under

352

operational conditions. The ensuing scaling effects would be further amplified in a future expansion of the

353

pilot BTES to full scale in the absence of a calibrated 3D model.

354

We further conclude that limitations of the analytical quasi-stationary BHE model utilized in FEFLOW

355

render simulations of short time charge and discharge pulses unreliable. For long term predictions, we

356

recommend the use of temporally aggregated flowrates and fluid temperatures in order to mitigate BTES

357

energy balance errors caused by simulation of short time BTES charge and discharge pulses.

358 359

Acknowledgements

360

the EUDP programme (project no. 64012-0007-1). The authors wish to thank Brædstrup Fjernvarme

361

(Brædstrup District Heating) for providing access to the operational data used in this work.

362 363

References

364

http://www.danskfjernvarme.dk/~/media/danskfjernvarme/videnom/publikationer/faktaark/faktaark_heat_ge

365

neration_in_denmark.pdf

366

[2] N. Giordano, C. Comina, G. Mandrone and A. Cagni, 2016. Borehole thermal energy storage (BTES).

367

First results from the injection phase of a living lab in Torino (NW Italy), Renewable Energy, 86, pp. 993-

368

1008.

369

[3] M. N. Fisch, M. Guigas, J. O. Dalenbäck, 1998. A review of large-scale solar heating systems in Europe,

370

Sol. Energy, 63 (6), pp. 355–366.

371

[4] M. Lundh, J.O. Dalenbäck, 2008. Swedish solar heated residential area with seasonal storage in rock:

372

initial evaluation, Renewable Energy, 33, pp. 703–711.

373

[5] T. Schmidt, D. Mangold, H. Müller-Steinhagen, 2004. Central solar heating plants with seasonal storage

374

in Germany, Sol. Energy, 76, pp. 165–174.

RI PT

350

M AN U

SC

The pilot BTES was constructed with support from the ForskEL programme (project no. 2010-1-10498) and

AC C

EP

TE D

[1] Danish District Heating Association, 2016. Heat generation in Denmark,

375 376

[6] H.-J.G. Diersch, D. Bauer, W. Heidemann, W. Rühaak, P. Schätzl, 2011. Finite element modeling of

377

borehole heat exchanger systems: Part 2. Numerical simulation, Volume 37, Issue 8, August 2011, pp. 1136–

378

1147.

379

[7] J. Nussbicker-Lux, 2012. The BTES project in Crailsheim (Germany) – Monitoring results, 12th

380

International Conference on Energy Storage - Innostock. 18

ACCEPTED MANUSCRIPT [8] B. Sibbitt, D. McClenahan, R. Djebbar, J. Thornton, B. Wong, J. Carriere, J. Kokko, 2012. The

382

performance of a high solar fraction seasonal storage district heating system – five years of operation,

383

Energy Procedia, 30, pp. 856–865.

384

[9] H. Bjørn, 2013. Borehole thermal energy storage in combination with district heating, Proc. European

385

Geothermal Congress 2013, June 3–7, Pisa, Italy, SG5-01.

386

[10] S. Hwang, R. Ooka & Y. Nam, 2010. Evaluation of estimation method of ground properties for the

387

ground

388

http://doi.org/10.1016/j.renene.2010.01.028.

389

[11] V. Wagner, P. Bayer, M. Kübert & P. Blum, 2012. Numerical sensitivity study of thermal response

390

tests. Renewable Energy, 41, 245–253, http://doi.org/10.1016/j.renene.2011.11.001.

391

[12] A. Casasso & R. Sethi 2014. Efficiency of closed loop geothermal heat pumps: A sensitivity analysis.

392

Renewable Energy, 62, 737–746, http://doi.org/10.1016/j.renene.2013.08.019.

393

[13] National well database (JUPITER), Geological survey of Denmark and Greenland (www.geus.dk)

394

[14] H.-J. G. Diersch, 2009. FEFLOW reference manual, WASY Gmbh, Institute for water resources

395

planning and systems research, Berlin.

396

[15] P. Eskilson, J. Claesson, 1988. Simulation model for thermally interacting heat extraction boreholes,

397

1988, Numer. Heat transf., 13(2), pp. 149-165.

398

[16] Balling, N., 2013. The Lithosphere beneath Northern Europe: Structure and Evolution over Three

399

Billion Years. Contributions from Geophysical Studies. Doctoral dissertation, Aarhus University. ISBN

400

9788790400620

401

[17] Doherty, J. 2010. PEST. Model-Independent Parameter Estimation. User Manual: 5th Edition,

402

Watermark Numerical computing, http://www.pesthomepage.org/getfiles.php?file=pestman.pdf.

403

[18] VDI, 2004. VDI 4640 Thermal use of the underground, VDI-Gessellschaft Energie und Umwelt (GEU),

404

Berlin.

heat

pump

system.

Renewable

Energy,

35(9),

2123–2130,

405

19

AC C

EP

TE D

M AN U

SC

source

RI PT

381

ACCEPTED MANUSCRIPT

RI PT SC M AN U TE D

•

EP

•

A finite element model of a pilot borehole thermal energy storage is developed. We determine improved thermal conductivity estimates using inverse modelling of distributed storage temperature measurements. Scaling effects are examined by comparing the calibrated thermal conductivity distribution to a simple thermal response test based estimate. The calibrated model yields improved storage energy balance predictions.

AC C

• •