Thermal effect of a borehole thermal energy store on the subsurface

Abstract

Background

The thermal effect on the subsurface of a large borehole thermal energy store (BTES) has been investigated by coupling measured rock properties with an enhanced FEFLOW simulation.

Methods

The finite element model has been validated against measured data from a 2-year operation period. The thermal changes in the subsurface have been predicted by simulation for a 30-year operation period. The model is based on three 80-m core sections drilled in Triassic carbonates, which have been analyzed in detail with respect to lithology, facies, and thermal and hydraulic parameters.

Results

The model shows thermal effects of the BTES on the subsurface at a distance of approximately 350 m after 10 years and a distance of approximately 850 m after 30 years of operation. At a distance of 100 m, the temperature of the subsurface rises by 2 K after 30 years.

Conclusions

The simulation describes the real BTES in an accurate manner and is suited for predicting the thermal changes in the subsurface for long-term operational durations.

Background

The present work intends to predict the thermal effect of a borehole thermal energy store (BTES) on the subsurface in the case of flowing groundwater by simulation with the software FEFLOW. Predicting thermal effects on the subsurface or grouting using finite element modeling has been done by other authors before (Wołoszyn and Gołaś [2013]; Rees and He [2013]; Wagner et al. [2013]; Zhang et al. [2011]), but the properties of the surrounding geological formation itself are generally simplified or focus on the grouting material only. Predicting thermal effects on the subsurface is essential for construction approvals by authorities and should not exceed a threshold value of 6 K ([VDI-Richtlinien 2011]). In this study, a detailed examination of the geology has been included into the FEFLOW model. The model has been further enhanced by integrating a newly developed add-on, which significantly reduces computation times and allows the coupling of TRANSYS with FEFLOW (Bauer et al. [2011]). The resulting model has then been validated against measured data spanning a 2-year period of operation of the BTES in Crailsheim.

Borehole thermal energy store in Crailsheim

A BTES is an underground heat exchanger system used to store large quantities of excess heat. A common application is the storage of heat from a combined heat and power plant or the capture of solar heat in summer for later use in winter. BTES consists of multiple arranged and interconnected borehole heat exchangers (BHEs) installed with a spacing of 2 to 5 m (VDI-Richtlinien [2001]). The thermal interaction of the individual borehole heat exchangers is essential for an efficient storage process.

The investigated BTES was constructed in 2008 in Crailsheim, Germany (Figure 1) (Riegger [2008]) and forms a major part of a solar district heating system (Bauer et al. [2008]). The BTES consists of 80 double U-tube BHEs, drilled to a depth of 55 m and separated at a distance of 3 m to each center. The well pad occupies a circular area of 30-m diameter. Each BHE is installed in a 0.13-m borehole and is cemented in place with thermally enhanced grout. The BTES is covered with 0.4 m of a thermal insulation material and 1.6 m of elevated soil from a nearby earthwork. The BTES in Crailsheim operates with pure water as working fluid. The BHEs are interconnected by pairs, with one BHE being situated more inside of the BHE grid and one more outside (Figure 2). The water flow direction is from the inner BHE to the outer during heat injection and the other way round during heat extraction.

Figure 1

Map showing the study area near Crailsheim.

Figure 2

Layout of Crailsheim BTES with position of groundwater monitoring wells and temperature measuring lances.

The BTES was monitored in detail during a 2-year operation period. The flow rate within the pipes from and to the heat storage, as well as the in- and outflow temperatures, is logged continuously during the operation phase of the system. From this data, the loading and the unloading heat amount can be calculated. Furthermore, 24 temperature sensors have been installed at 12 different depths on measuring lances M1 and M31 (see Figure 2). Both measuring lances record temperatures of the solid underground up to a depth of 80 m.

In addition to monitoring the flow rates and fluid temperatures, groundwater temperatures and levels are measured in several groundwater monitoring wells. GWM2, GWM3, and GWM4 are 80 m below ground surface (aquifer in the Upper Muschelkalk), and in GWM2a, GWM3a, and GWM4a, wells are 20 m below ground surface (aquitards in the Lower Keuper) (see Figure 3).

Figure 3

BTES Crailsheim well section - stratigraphy, facies, and hydrogeology parameters. (a) Clay/dolomite. (b) Peloidal grainstone. (c) Mud/wackestone. (d) Bioclast grainstone.

Geological setting

The study area is part of the Western Germanic Basin. During Triassic times, the Germanic Basin was situated in a peripheral position towards the Tethys spreading center and also affected by syndepositional tectonism. Generally, two main systems governed the basin's tectonic evolution. The first, mainly in the eastern and southern parts of the basin, was controlled by reactivated Variscan structures; the second developed in the NW part of the basin related to the North Sea rifting belt (Feist-Burkhardt et al. [2008]). During Middle Triassic times, the Variscan structures were reactivated and thereby influenced the sedimentary processes and tectonic subsidence within the study area. Today's fracture system resulted from the stress field that has developed since the Late Cretaceous. Main fracture directions are NNE-SSW and NE-SW oriented. These two main directions were also measured in the quarries Heldenmühle and Neidenfels near Crailsheim showing a dominant NE-SW orientation.

The study site is located in the Hohenloher Ebene which composes the largest Muschelkalk karst area in Germany as part of the South German Scarplands. The Hohenlohe area is dominated by sedimentary rocks of Triassic age. The stratigraphic formations are (from the oldest to the youngest): the Middle Triassic Muschelkalk subdivided into Lower Muschelkalk (Wellenkalk), Middle Muschelkalk, and Upper Muschelkalk (Hauptmuschelkalk), and the Upper Triassic Keuper which is subdivided into Lower Keuper (Lettenkeuper), Middle Keuper (Gipskeuper), and Upper Keuper. Bedding is nearly horizontal with a dip of about 3° SE.

The Muschelkalk in the Hohenloher Ebene reaches a total thickness of 250 m with Lower, Middle, and Upper Muschelkalk series of about 80-m each. The Upper Muschelkalk is overlain by a 10- to 30-m-thick series of sandstones and claystones of the Lower Keuper.

Hydrogeology

Hydrogeologically, drill logs and flow meter measurements performed in all groundwater monitoring wells suggest multiple aquitards separated by layers of clay in the topmost 50 m of the groundwater monitoring wells. Fracture-controlled aquifers in the drilled zone exist in depth deeper than 50 m. Groundwater level measurement data indicates a flow direction of NE-SW for the fracture-dominated aquifer. This direction is supported by the general NE-SW strike of the Upper Muschelkalk. The flow direction of the aquitards can only be inferred on a regional scale to be SSW-NNE.

No regional groundwater flow rates were available for this area. Flow meter measurements in the wells could not detect any significant flow rates, neither through fractures nor through matrix. However, studies by the Bavarian Environment Agency ([2007]) in southwest Germany suggest hydraulic conductivities of 1.0 × 10⁻⁵ to 1.0 × 10⁻⁴ m s⁻¹ for the Middle and Upper Muschelkalk and 1.0 × 10⁻⁶ to 1.0 × 10⁻⁴ m s⁻¹ for the Lower Keuper. Pumping and recharge tests conducted at the six monitoring wells returned transmissivities in a range of 1.9 × 10⁻⁷ m² s⁻¹ to 5.1 × 10⁻⁵ m² s⁻¹. Storage coefficients (S_S) returned from hydraulic tests in all groundwater monitoring wells imply unconfined groundwater conditions (S_S 8.0 × 10⁻³ m⁻¹ to 2.0 × 10⁻² m⁻¹) for the aquifer in the Upper Muschelkalk and confined groundwater conditions (S_S 6.0 × 10⁻⁵ m⁻¹ to 3.0 × 10⁻⁴ m⁻¹) for the aquitards in the Lower Keuper.

Well logs

The geological formations penetrated by the 55-m-deep BTES belong to the Upper Muschelkalk and the Lower and Middle Keuper. The Middle Muschelkalk was reached by an adjacent exploration well at about 100-m depth. Groundwater monitoring wells GWM2, GWM3, and GWM4 have been cored over the whole length (80 m each); thus, the lithology of all drill cores can be logged in detail (Figure 4). Stratigraphic correlation of all three groundwater monitoring wells is very good, whereas most marker horizons commonly differ in less than 2 m in intersected depths. Additionally, the petrography of GWM3 has been examined in detail using thin sections and acetate peels of representative stratigraphic intervals. Stratigraphic correlation of all three wells was possible since they are located within an area of about 150 m² in horizontally layered formations.

Figure 4

Borehole heat exchanger modeled as 1D line element in a 3D finite element mesh.

GWM3 starts with coquina beds of the Crailsheim Beds (Trochitenkalk Formation, mo1) in 80-m depth, extending up to about 72-m depth below ground surface. The Crailsheim Beds are characterized by a limestone-marl alternation with peloidal packstones. The following Bauland Beds range up to about 59 m and consist of up to about 63 m of clayey limestone with intercalated microcrystalline packstones. The topmost 4 m of the Bauland Beds comprises bioclastic packstones, grainstones, and rudstones.

The following Meißner Beds (mo2) form the lowest subunit of the Meißner Formation ranging up to about 47 m. This unit primary consists of packstones, rudstones, and coquina. Clay horizon β, a regional marker horizon, is located at about 53-m depth. The Meißner Beds are separated from the overlaying Künzelsau Beds (mo3) by the clay horizon γ, a second marker horizon in this area. Künzelsau Beds consist of micritic limestones with intercalated microcrystalline packstones and gypsum nodules. A third marker horizon, clay horizon δ, is located at about 44-m depth. Fine-grained coquina beds in 41-m depth are followed by bioclastic and peloidal grainstones up to about 35-m depth. The topmost 10 m of the Meißner Formation are formed by limestone with intercalated dolomitic marls, micritic limestone, and coquina. The boundary between the Upper Muschelkalk and the Lower Keuper is located at about 24-m depth. The Lower Keuper starts with dolomitic marl, sandstone, limestone, and intercalated claystone. At about 15-m depth, massive sandstones are present, followed by layers of marl and dolomite. The topmost 4 m of the well is formed by anthropogeneous clay.

Modeling

In general, the geological and hydrogeological settings at the site are very suitable for validating the recently developed FEFLOW add-on in conjunction with a BTES (Bauer et al. [2011]). Geological and hydrogeological settings such as the low, but still measurable, groundwater flow through the BTES and the aquifer below the BTES with its higher flow rate enable the investigation of the thermal effect of a BTES on the subsurface and the verification of the developed simulation software.

The developed simulation software is an add-on to the finite element flow and transport modeling tool FEFLOW (Zheng [2007]; Diersch et al. [2011a, b]). The add-on is capable of taking groundwater flow into account for large borehole heat exchanger fields. Borehole heat exchangers are modeled in FEFLOW as one-dimensional line elements in a three-dimensional finite element model of the underground. Heat and mass transport processes inside the BHE are modeled based on thermal resistance and capacity models (Bauer et al. [2011]). Discretizing the surrounding ground with a relatively coarse mesh (compared to a three-dimensional discretization of a BHE) means that the number of finite elements can radically be reduced; this saves computation time during simulation. BTES can be modeled and simulated in an efficient and variable manner with this technique. Figure 4 shows a one-dimensional line element representing a borehole heat exchanger in the three-dimensional finite element modeling of the surrounding ground.

Methods

The thermal conductivity and thermal diffusivity, permeability, porosity, and density have been measured on 76 representative samples of drill cores from GWM3. Conventional methods for determining permeability commonly return an integrative result that refers to the whole sample. A gas pressure permeameter permits a numerous number of pointwise measurements, which can be assigned to the individual points of a rock sample.

Laboratory measurements

Permeability measurements were performed with a combined columnar and minipermeameter after Jaritz ([1999]) and Filomena et al. ([2013]). The minipermeameter is a compressed air-driven gas permeameter. It allows a fast and precise measurement of the permeability at the surface of a rock sample in comparison to other permeability-measuring techniques. The sample's surface can be plane or curved (e.g., a drill core), but has to be smooth. The measurement is fast and nondestructive. Multiple measurements have been done on each sample (scattered over all suitable surfaces), and the readings have been averaged to get one value per sample.

Furthermore, a columnar permeameter is used to validate the results of the minipermeameter. The columnar permeameter requires test slices that are sawed off from selected drill cores. It is necessary to consider moisture and temperature of the ambient air when measuring permeability. For this reason, the measuring accuracy is significantly improved by using ovendry samples. Similarly, earth-moist samples cannot be measured. A detailed study on the technique and accuracy of gas permeameters, including the device used for this study, can be obtained from Filomena et al. ([2013]).

Porosity and density have been measured using a commercial gas-driven pycnometer combined with a displacement-measuring technique that utilizes a fine-grained powder as quasi-fluid for the determination of the sample volume. The method allows fast but highly accurate measurements (with 0.03% accuracy and 0.02% reproducibility) of high quantities of samples, but only returns the gas-effective porosity.

Thermal conductivity and thermal diffusivity were measured on ovendry samples applying the optical scanning method after Popov et al. ([1999]). The method allows relatively fast measurements of samples showing differing quality with an accuracy of 3% (Bär et al. [2011]).

Model setup

The heat and mass transport model was created using the numerical simulation modeling software for subsurface flow and transport processes FEFLOW 6.0. The model dimension extends over 1.5 km in the downstream direction and 500 m in the upstream direction along the aquifer. The width of the model was set to 1 km, with the BTES being placed in the middle. The depth of the model was set to 150 m with the BTES ranging from 0 to 55 m in depth.

Model layers and their thickness have been chosen individually to fit the geological realities. However, as some sections have an intensely interbedded structure with layers often changing every few centimeters (e.g., the Bauland Beds, Meißner Beds, and Künzelsau Beds), they could not be modeled in detail. Consequently, multiple single geological layers had to be merged into reasonable model layers. The maximum discretization level for the vertical element heights was set to 1 m and the minimum to 10 m. Generally, layers in the vicinity of the BTES (i.e., the topmost 75 m) have a thickness of 1 to 4 m while layers far below the BTES have a thickness of 5 to 10 m. The upper part (about 5 m) of the boreholes has a different diameter and is filled with a different grouting material (Riegger [2008]). This was not modeled as only a minor effect on the thermal behavior of the store is expected. The covering of the BTES, consisting of a thermal insulation material (0.4 m) and elevated soil from a nearby earthwork (1.6 m), was modeled by two additional continuous layers that taper off to 0.01 m each outside the BTES. The predefined mesh size changes with distance towards the BTES. Element sizes around each individual BHE are finest with horizontal triangle areas of 0.02 to 0.08 m², increasing to 0.1 to 0.4 m² in the interspace between each BHE and increasing further to 0.5 to 2 m² in an oval-shaped area that encloses the BTES and points in the downstream direction. The model has its coarsest elements (200 to 500 m²) in the remaining area. The total model consists of 57 layers and about 1,445,000 finite elements (Figure 5).

Figure 5

Model meshing and layer dimensions of the 3D finite element mesh. Vertical cut through the middle of the thermal storage along the aquifer flow direction; vertical exaggeration, ×5.

Groundwater flow directions were defined by setting suitable hydraulic heads as a boundary condition of the first kind at the model borders at the depth of the aquitards and the aquifer. The hydraulic gradient of the aquitards was set to 0.01, while the aquifer has a value of 0.0045. These values were given by the Regional Authorities for Geology, Raw Materials and Mining as typical values for the site.

Based on the available measuring results, hydraulic conductivities were assigned to the individual model layers. Measurements were averaged when multiple measurements exist within the range of an individual model layer. For layers deeper than 80 m, data from unpublished internal studies were used.

After setting up the hydraulic parameters, the initial hydraulic conditions for the model were simulated in steady flow mode (Table 1).

Table 1 Layer dimensions and input data of the FEFLOW model

Subsequently, the unsteady heat transport simulation was then set up by assigning porosity, thermal conductivity, and heat capacity to each layer of the flow model. Measurement data of the thermal conductivity and heat capacity of the solids were measured on dry samples and thus have to be converted to fully saturated conditions. In FEFLOW, this can be achieved by inputting the thermal properties of the fluid (i.e., the groundwater), which were set to 4.19 × 10⁺⁶ J m⁻² K⁻¹ for the volumetric heat capacity and 0.6 W m⁻¹ K⁻¹ for the thermal conductivity.

Porosity measurements using a gas pycnometer return the gas-effective porosity and thus have to be converted to the water-effective porosity first. As no water-based measurements where possible due to decomposition and/or swelling of the samples, known published porosity data from similar rocks (Götz and Lenhardt [2011]) were used as a reference. By decreasing all porosity measuring results by a factor of 5, the resulting porosity magnitude accords quite well with the published values, while the porosity varieties of the measured samples are preserved.

The initial temperature distribution of the model has been calculated by using a fixed surface temperature of 10.5°C (Blocon [2014]) and a heat inflow of 0.065 W m⁻² (Pollack et al. [1993]) at the base of the model. The resulting temperature gradient is about 3.5 K per 100-m depth.

The 80 BHEs have been integrated as a boundary condition of the fourth kind, which is capable of simulating injection and extraction at nodal points. Grout and pipe properties have been set according to data provided by the manufacturers. After setting up all required BHE parameters, all of the BHEs were interconnected pairwise using FEFLOW's internal BHE Integrator.

At first, the transient simulation was conducted over 2 years in order to validate FEFLOW against measured data from a 2-year period of operation. The BTES was only charged during the first 2 years. Then the measured BHE thermal load was set as a boundary condition to the model (see Figure 6). Then the computed and measured temperatures at measuring lances M1 and M31, as well as groundwater temperatures at GWM2a to GWM4a and GWM2 to GWM4, were compared.

Figure 6

Measured BHE thermal load set as a boundary condition to the model for validation purpose.

Finally, the thermal changes in the subsurface were simulated over a time span of 30 years. The semiannually alternating operation of the BTES was simulated according to the planning with inflow temperatures of 75°C in the summer months (April to September) and 20°C in the winter months (October to March). Simulation duration using a common workstation (Quadcore with 2.66-GHz clock speed and 24-gigabyte internal memory) took about 12 h.

Results and discussion

Figure 7a,b depicts the comparison between measured and computed underground temperatures at measuring lances M1 (middle of the BTES) and M31 (margin of the BTES) of the 2-year validation run. Figure 7c,d depicts the comparison between measured and computed groundwater temperatures at groundwater monitoring wells GWM2a to GWM4a (20-m depth) and GWM2 to GWM4 (80-m depth).

Figure 7

Measured and computed underground and groundwater temperatures. Measured and computed underground temperatures at measuring lances M1 (a) and M31 (b) and measured and computed groundwater temperatures at wells GWM2/a, GWM3/a, and GWM4/a at 20-m depth (c) and 80-m depth (d).

It appears as though both, the temperatures of the solid ground and the groundwater, can be properly computed with the FEFLOW model. Minor differences in results should be accepted because the measured BHE thermal load set as a boundary condition to the model is slightly too high. This is because it includes the unknown thermal losses of the 300-m-long pipes to and from the BTES causing underground temperatures to be computed higher than normal at the end of the second year. Furthermore, it is not possible to install the BHE and the measuring lances absolutely vertical. Hence, an uncertainty of measurement due to inexact sensor positioning must be considered. The measured groundwater temperatures at GWM2a and GWM3a show an unexpected course in March 2010 and the end of June 2010, respectively. This may be explained by the introduction of surface water which led to an erroneous reading.

Altogether, it can be stated that the simulation model describes the real BTES in an accurate manner. It is suited for predicting the thermal changes in the subsurface for an operational duration of 30 years.

Figure 8 shows the computed temperature profiles in the subsurface in vertical cross sections parallel to the aquifer flow direction after 3, 10, and 30 years of operation. The results of the simulation show a clear impact of the BTES on the aquifer, while the aquitards are only slightly affected. Apparently, the transmissivity of the aquitards is too low to significantly impact the temperature profile.

Figure 8

Computed temperature profiles in the subsurface in vertical cross sections parallel to the aquifer flow direction. After 3, 10, and 30 years of operation; vertical exaggeration, ×5.

After 3 years of operation, the advective-caused heat transport in the main flow direction is almost invisible. However, after 10 years, an obvious temperature propagation up from 350-m distance to the BTES is visible in the aquifer. The propagation increases up to 850 m after 30 years of operation. Nevertheless, the temperature change in this distance is below 1 K. In 100-m distance to the BTES, in the aquifer flow direction, the computed temperature rise at the most affected depth (71 to 76 m) is 2.1 K.

The temperature profile inside the BTES is only slightly affected by advective effects (Figure 9). This implies a low impact of the groundwater flow on the heat storage efficiency. Furthermore, results show a distinctly higher temperature spreading in the deeper parts of the BTES between a completely charged and discharged condition caused by the higher thermal conductivities of the Bauland Beds. By volume, this allows increased heat storage, albeit the heat loss increases, too.

Figure 9

Computed temperature profiles of the BTES in the 30th year of operation. Completely charged (a) and completely discharged (b); left: vertical cross section parallel to the aquifer flow direction; right: horizontal cross section at approximately 40-m depth; temperature scale in degrees Celsius.

Conclusions

Three-dimensional modeling and simulation of a BTES located predominantly in Middle Triassic carbonates with multiple aquitards at the upper part of the BTES and one aquifer below the BTES were carried out in order to determine the thermal effect of the BTES on the subsurface. As the BTES was modeled in detail with real technical and measured geological and hydrogeological data, the resulting model is of high quality. This was confirmed by validating the computed transient temperature changes of the solid underground and the groundwater against measured data from a 2-year period of operation.

Possible thermal changes in the surrounding subsurface have been predicted by simulation for an operation duration of 30 years. The model shows thermal effects of the BTES on the subsurface at a distance of approximately 350 m after 10 years and a distance of approximately 850 m after 30 years of operation. At a distance of 100 m, the temperature of the subsurface rises by 2 K.

All considered, the examined simulation program, FEFLOW, is a powerful tool for projecting the thermal effects of seasonal underground heat storage systems on the subsurface in a time-efficient manner. Thus, FEFLOW is of importance for authorities and investors for granting and planning similar systems in future projects.