The wider SKEWS project area is penetrated by several geotechnical drillings and multiple groundwater monitoring wells. This data was extracted and evaluated from the Hessian drilling database (HLNUG 2020). Hence, a first-order geologic and hydrogeologic model concept was developed and implemented into FEFLOW.
Geologic and hydrogeologic exploration
A total of 325 drilling datasets from the archive dating from 1950 to 2020 were provided by the Hessian agency for Nature Conservation, Environment and Geology (HLNUG) and evaluated to gain additional information about the local distribution of the different lithologies (Fig. 3a)
The evaluation of the drilling database allowed more precise localization of the suspected fault zone. After surveying the litho-logs, the fault zone was assumed in the eastern part of Lichtwiese Campus close to the planned MD-BTES. The information gained from the drilling data was imported into the suite ESRI ArcGIS. An interpolator was then used to create a raster file representing each geologic unit’s depth map. Since only shallow drillings were available for the sedimentary Rotliegend in the east, a thickness of 40 m on top of the crystalline basement was assumed based on a minimal thickness of Rotliegend Sediments observed in drillings close to the suspected fault 350 m north of Lichtwiese campus.
To parametrize the groundwater model, data of 42 monitoring wells were used to generate transient contours with a monthly resolution (Fig. 4). The data were evaluated from 2000 to 2020. A monthly average of the hydraulic head was calculated for each monitoring well to integrate the dynamics of the hydraulic head.
Average groundwater contour lines for each month for the period 2000–2020. In red the FEFLOW model boundary, derived from the groundwater flow direction
The resulting single aquifer piezometric map shows a south–north direction of groundwater flow, which deviates to the west at the North of Campus Lichtwiese towards the URG (Fig. 4). The depth of the phreatic surface displays an average of approx. 5 m below surface. The model boundary for the FEFLOW model was based on the groundwater contours and the flow lines of the groundwater table. The goal of the chosen model boundary is to minimize the influences of the model boundary conditions on the BTES by providing sufficient separation between the boundary conditions and the MD-BTES. Additionally, the model was extended to the North to investigate the potential thermal influence on nearby water bodies like the Woog pond (Fig. 4).
Finite element model
The construction of the SKEWS system is planned with up to 4750-m-deep pilot borehole heat exchangers. This study was finished in the period of the drilling operation. Therefore, a simplification of the complicated geological conditions that were entcountered had to be made. A preliminary study on the storage test operation performed within the SKEWS project revealed a negligible thermal impact on the aquifer. However, a step-wise enlargement to 37 BHE in total is envisaged after the successful finalization of the SKEWS project following the pre-calculated storage models. Consequently, the investigations within this study consider the regular seasonal operation of a storage system that comprised all 37 BHE. The BHE array has a hexagonal shape with a BHE spacing of 5 m (Fig. 5). Due to their enhanced robustness in higher depth, coaxial pipe systems are used in SKEWS. In theory, coaxial BHE is advantageous because the flow direction can be reversed to reduce heat loss during loading and unloading cycles (Welsch et al. 2016), which could not be implemented due to the limitations of the FEFLOW-internal BHE solution. Furthermore, a thermal insulation of the upper part of the BHE, for example with the tool BASIMO (Schulte et al. 2016b, c), is not implemented in the model for reasons of simplification. The parametrization of each BHE is displayed in Table 1.
Location and layout of the planned MD-BTES with the three building stages and additional groundwater monitoring wells
The FEFLOW model covers an area of 2 km2. Its depth was set to 1000 m to prevent an influence of the lower boundary conditions on the BHE. A structured mesh of prismatic finite elements was used (Fig. 6), consisting of 87 mesh layers. Moreover, the mesh was refined around the BHE positions and closer to the fault zone to obtain better calculation results, where high-temperature variations or steps in the material properties due to the fault displacement are likely to occur (Fig. 6). The BHE are implemented as a particular type of boundary condition, at which it is possible to add or extract a defined amount of heat into the model (Diersch 2014). A dual continuum approach was applied for BHE simulations. It features solving of fluid and heat transport processes in the natural subsurface with a 3D finite element approach. The BHE, in contrast, are computed using an analytical solution (Diersch et al. 2011b; Bauer 2011; Eskilson and Claesson 1988), which is coupled to the finite element mesh along with 1D finite element representations. To ensure an optimal coupling of the BHE solution to the 3D model, the nodes surrounding the BHE (help nodes in Fig. 6) were placed at an optimal node distance (Diersch et al. 2011b; Nillert 1976). The approach was chosen due to its robustness and reasonable accuracy for long-term simulations (Diersch et al. 2011a, b). In a previous study, the simulations agreed with measured data (Mielke et al. 2014).
Model domain with the integrated BHE and the discretization scheme around every BHE node
The geological drilling data were imported from the layer contour map as a point set and interpolated into a meshed surface in FEFLOW. In total, five horizons were imported: the ground surface from a digital elevation model, the base quarternary, the base Rotliegend and the base granodiorite weathering zone. The displacement was assumed to be constant with 30 m for the whole fault. Due to the postulated paleozoic weathering of the crystalline, the weathering zone of the granodiorite was continued underneath the Rotliegend, but with reduced thickness (Fig. 7).
FEFLOW 3D model with the distribution of the lithologies
This geological model was created using the generated surfaces as a unit basis. One disadvantage of the layered approach is that layers that the fault would cut off still propagate through the whole model domain. The resulting errors were minimized by fitting the respective material properties of the corresponding correct lithology to these elements.
Hydraulic properties
The model is considered to be fully saturated. Three boundary conditions (BC) were assigned. Firstly, hydraulic head BCs (first-order Dirichlet BC) were applied to the northern and southern boundaries. Values for the hydraulic head BC were derived from the groundwater head data as a mean of the values from Fig. 4. Groundwater recharge of 170 mm a−1 was assumed for the top of the model, considering an approximate degree of artificial soil sealing of 10% at the Lichtwiese campus (Beier 2008). Hydraulic conductivity parameters (Table 2) for the sedimentary units and the granodiorite weathering zone were extracted from pumping tests at the project site. The hydraulic conductivity for the basement shown in Table 2 is based on formation permeability data for crystalline rocks (Stober and Bucher 2007).
Thermal properties
For the temperature, two first-order Dirichlet Boundary conditions were applied to the top and bottom of the model. Accordingly, the temperature at 1000 m depth was assumed to be 50 °C (Rühaak et al. 2012), corresponding to a geothermal gradient of 4 °C per 100 m. For the surface temperature, the average temperature in Darmstadt for the years 1981 to 2010 of 10.1 °C was used (DWD 2019). Thermophysical rock properties for the different geologic units (Table 2) were derived from outcrop analog studies of regional outcrops (Bär et al. 2019; Weinert et al. 2020).
Thermal conductivity and volumetric heat capacity are temperature-dependent parameters. The thermal conductivity decreases with increasing temperature, whereas the heat capacity shows the opposite behavior (Mottaghy et al. 2008; Schatz and Simmons 1972; Ahrens 1995). Hence, both parameters were temperature corrected prior to simulation (Vosteen and Schellschmidt 2003; Bär 2012).
Parameter fitting
A parameter fitting for the hydraulic conductivity was conducted, to adjust the hydraulic properties of the model units so that the simulated hydraulic heads fit the field data. The parameter fitting was done with Pest (Doherty 1994) using a pilot point approach (Doherty 2003; Fienen et al. 2009).
Implementation of the permeable fault
With reference to the pre-drill geological and geophysical knowledge, the hydraulic fault properties were initially set to low permeability. While drilling it became obvious that the fault provides a reasonable hydraulic conductivity (e.g., large circulation losses). Therefore, a broader range of fault hydraulic properties was implemented. Exact results for fault permeability can be expected after evaluation of the drilling and pumping data.
The interpreted fault was implemented into the model using a discrete feature element. Discrete feature elements are 1D or 2D objects added to an existing 3D FEFLOW mesh to simulate features with specific flow properties like fractures or faults (DHI WASY 2009). Here, a discrete feature element applying Darcy’s law for the fault internal flow simulation was used. As concluded from the evaluated drilling data, the implemented fault zone only interferes with one-half of the borehole storage system (Fig. 8)
Position of the fault in the model in the 2D view and 3D view
Variation of fault parameters
Since the fault zone characteristics are not well known, a parameter variation study was conducted to investigate the effects of different hydraulic conductivities and thicknesses of the fault zone on the temperature field and the storage system performance. Accordingly, the respective properties of the discrete feature elements were varied. In total, 15 models with fault zone thicknesses of 0.5 m, 1 m, 2 m, 3 m, and 4 m and hydraulic conductivities ranging from 10−3 m s−1and 10−4 m s−1 to 10−5 m s−1 were simulated. Fault thicknesses were chosen to get a reasonable contrast between the simulation results.
Borehole heat exchanger operation scenario
To obtain comparable results, a simplified loading and unloading scheme was applied for every simulation (Fig. 9). During charging periods, the inlet temperature was set to 90 °C. This temperature can be supplied by solar thermal collectors and describes the upper-temperature limit of PE-X pipes (Welsch 2019). During the heat extraction period, the inlet temperature was set to 30 °C to ensure the supply of low-temperature heating systems. Each loading and unloading cycle ran for half a year (183.5 days). The total simulation time is 30 years (10950 days). BHE were connected in parallel to provide a constant inlet temperature for each of them. For the whole simulation time, the flow rate for each BHE was set to 2 l s−1. According to some preliminary considerations, this value represents a good compromise between a comparably low pressure drop in the BHE and a high heat exchange rate.
Exemplary inlet and outlet temperatures for the whole BHE array
Data analysis
Multiple performance indicators were calculated to assess the efficiency of the different models. Firstly, the mean outlet temperature was derived for the whole BHE array as an average of all BHE. Heat is transferred continuously due to the temperature difference between the heat transfer fluid and the surrounding rock. The heat rate (dot{Q}), which is the heat exchanged between the heat carrier fluid and the surrounding rock, was calculated with Eq. 1:
$$dot{Q}=Delta Tcdot {(rho c)}_{f}cdot dot{V,}$$
(1)
where ∆T is the difference between inlet and outlet temperature of the heat carrier fluid with a flow rate of (dot{V}) through the BHE array. (ρc)f is the volumetric heat capacity of the fluid. The fluid parameters were assumed to be constant for the whole temperature range to simplify the system. The total heat stored or extracted could be calculated by integrating exchanged heat over a charging or discharging cycle (1 year). The ratio between the amount of heat stored and extracted is defined in Eq. 2 as the storage efficiency:
$$mu =left|frac{{Q}_{E}}{{Q}_{s}}right|,$$
(2)
with QE being the total amount of heat extracted and Qs the total amount of heat stored during a period of 1 year.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Disclaimer:
This article is autogenerated using RSS feeds and has not been created or edited by OA JF.
Click here for Source link (https://www.springeropen.com/)
