Assessment of a medium-deep borehole thermal energy storage site in the crystalline basement: A case study of the demo site Lichtwiese Campus, Darmstadt

The use of intermittent renewable energy sources is of great importance for reducing greenhouse gas emissions. Medium-deep borehole thermal energy storage systems (MD-BTES) represent an economic solution. At the Technical University of Darmstadt, Germany, an MD-BTES consisting of three 750 m deep borehole heat ex-changers was constructed as a demonstrator. Before construction, a comprehensive dataset consisting of electrical conductivity tomography profiles, gravity measurements, 2D seismic profiles


Introduction
Challenges such as the climate crisis and the recent increases in energy costs underline the importance of developing an efficient heat supply scheme.The European Green Deal targets a reduction in greenhouse gas emissions of 55% compared to 1990 levels until 2030 and climate neutrality by 2050 (European Commission 2019).As the building sector accounts for up to 40% of total European energy consumption and, thus, 36% of CO 2 emissions (Economidou et al., 2020), buildings are associated with a significant energy-saving potential (Tsemekidi-Tzeiranaki et al., 2020).
Besides increasing buildings' efficiency, e.g., by modern insulation, integrating renewable heat sources plays a key role in a renewable heat supply scheme.A central pillar of replacing fossil fuels in the heat supply is seasonal heat storage, which can buffer the often irregular or weatherdependent heat production, for example, from solar thermal energy or waste heat (Xu et al., 2014;Yang et al., 2021;Guelpa and Verda 2019).A promising approach for underground heat storage is borehole thermal energy storage (BTES) (Gehlin 2016;Lanahan and Tabares-Velasco 2017;Lizana et al., 2018).BTES utilizes the subsurface as a heat storage medium with conductive borehole heat exchangers (BHE).However, the already well-established shallow BTES (Fisch et al., 1998;Schmidt et al., 2004;Lundh and Dalenbäck 2008;Bauer et al., 2010;Sibbitt et al., 2012;Xu et al., 2014;Tordrup et al., 2017) have a large surface footprint and are often located in aquifers used for groundwater extraction (Bär et al., 2015;Welsch et al., 2016).Thus, they conflict with the water supply, as the increased temperature can affect negatively the groundwater (Griebler et al., 2015;DGG and DGGT 2015) For this reason, after extensive numerical studies (Welsch et al., 2016;Welsch et al., 2018;Schulte et al., 2016), a medium-deep BTES (MD-BTES) was constructed at the Technical University of Darmstadt as part of the SKEWS (Seasonal Crystalline Borehole Thermal Energy Storage) project.MD-BTES targets crystalline rocks with low hydraulic conductivity and high thermal conductivity at depths of up to 1000 m (Welsch et al., 2018).MD-BTES systems have a reduced surface space requirement compared to shallow BTES, making them especially suitable for urban areas and associated heating networks.The generally lower permeability at greater depth reduces convective heat transport by groundwater (Welsch et al., 2016).Furthermore, the near-surface area of the storage can be thermally isolated to ensure that the heating of the near-surface groundwater is minimized.Additionally, heat losses are reduced due to the higher temperature in deeper rock formations (Fig. 1).
As 33% of the rock exposed on the surface worldwide is constituted of crystalline rocks (Blatt and Jones 1975), there are many potential urban centers to expand this technology.To the knowledge of the authors, this project constitutes the world's first MD-BTES in crystalline formations, meaning that experience from previous projects cannot yet be drawn upon.
In general, the drilling phase accounts for the largest proportion of the investment costs of geothermal projects (Anderson and Petty 2006;Lukawski et al., 2014).Such costs and risks associated with drilling are already higher for MD-BTES compared to shallow BTES due to the increased drilling depth.The risk associated with drilling will increase dramatically in cases of insufficient geological knowledge.For MD-BTES, the local geological setting influences storage performance, and also the environmental impact, which is strongly dependent on rock-specific parameters (Rühaak et al., 2015).For instance, hydraulic conductivity and thermal conductivity (Welsch et al., 2016).Therefore, geologic uncertainties make it difficult to predict the BTES performance before drilling, complicating the proper dimensioning of all system components (Schelenz et al., 2017;Vienken et al., 2016).
Thus, geological uncertainty should be reduced before the drilling phase of BTES through geologic and geophysical exploration (Hesselbrandt et al., 2021).Careful surveying is particularly important for crystalline environments, often characterized by complex anisotropic and heterogeneous conditions (Bertrand et al., 2021;Bossennec et al., 2022).Previous publications on exploration in crystalline settings are mainly related to deep geothermal systems (Hloušek et al., 2015;Lüschen et al., 2015;Darnet et al., 2020;Rosberg and Erlström 2019;Frey et al., 2021), ore prospecting (Malehmir et al., 2012;Evrard et al., 2018) or engineering geology applications (Danielsen and Dahlin 2009).However, only very limited reports are available for exploring borehole thermal energy storage systems (Hesselbrandt et al., 2021) and no study has dealt with MD-BTES before.
This paper describes the procedure for exploring the MD-BTES site in Darmstadt.The multidisciplinary approach used is explained, emphasizing issues encountered and identifying possibilities for improvement during the construction of new MD-BTES systems.The geological situation based on the exploration results is explained and possible effects on the storage efficiency of the MD-BTES system are briefly outlined.

Geological setting
The site investigated in this study is located at the northern edge of the Odenwald crystalline complex, the largest outcrop of the Mid-German Crystalline High (Stein 2001), as well as at the north-eastern margin of the Upper Rhine Graben (URG).The crystalline rocks of the Odenwald are usually subdivided into two main petrogenetic units based on the individual intrusional and metamorphic history.These two main units of the western Bergsträßer Odenwald and the eastern Böllstein Odenwald, are separated by the Otzberg Shear Zone (Fig. 2a), a sinistral fault system (Krohe 1992;Stein 2001;Will et al., 2017) on a regional scale.The east of Darmstadt is part of the northernmost Bergsträßer Odenwald, dominated by large plutonic complexes intruding into a metamorphic volcano-clastic series (Altherr et al., 1999).The Bergsträßer Odenwald is divided into three subunits: the Frankenstein Complex (Unit I), the Flasergranitoid Zone (Unit II), and Fig. 1.Characteristics of MD-BTES in comparison to conventional shallow BTES.The MD-BTES targets the generally low permeable crystalline rock (in gray), which reduces negative effects on shallow groundwater aquifers (blue) in comparison to shallow BTES systems as the upper section of the storage can be thermally insulated (modified after Sass et al., 2016).the southern Odenwald (Unit III), resembling the Weschnitz Granodioritic Pluton and Tromm and Heidelberg Granitic Plutons (Krohe 1991).
The central Frankenstein Complex, which is relevant for the project site, (Unit I of the Bergsträßer Odenwald) consists predominantly of subduction-related gabbros, diorites, metasediments, and metabasalts.Radiometric dating of the gabbro revealed an intrusion age of 362 ± 7 Ma (Kirsch et al., 1988).Occurrences of granodiorites, diorites, and granites in the Frankenstein Complex are restricted to the Darmstadt area.These felsic to intermediate intrusions have K-Ar cooling ages of 340 ± 4 Ma and biotite cooling ages of 333 ± 4 Ma indicating a possible second magmatic cycle after the gabbro intrusions of the Frankenstein complex (Kreuzer and Harre 1975;Dörr and Stein 2019).South of Darmstadt, NE-SW striking metamorphic amphibolite and metapelite units are outcropping (Fig. 2a, b).The Dieburg Metagranite in the northern Frankenstein Complex was dated to 540 ± 8 Ma, recently (Stein et al., 2022).This study indicated for the first time the presence of Cadomian relics in the Variscan orogen.
The thickness of the Rotliegend deposits decreases from 250 m in the northern Sprendlinger Horst to zero in the south of Darmstadt (Kowalczyk 2001;Mezger et al., 2013).Overall the post-erosional sedimentary thickness is controlled by tectonics.The prevalent deposit of the Langen-Formation consists of fluvial sediments (Lippolt and Hess 1983) that are locally intersected or covered by basalts and andesites (Marell 1989;Fei et al., 2021).The Permian basalt is bound to SSW-NNE-striking structural elements of the Sprendlinger Horst (Fahlbusch 1980;Kowalczyk 2001).Lithostratigraphy is correlated with the volcanic deposits of the Donnersberg Formation in the Saar-Nahe Basin.They were dated to 290 ± 6 Ma (Lippolt and Hess 1983).These volcanic units often form the base of the Rotliegend lying directly on top of the crystalline basement.However, locally, sedimentary units can be found directly overlying the basement (Fei et al., 2021;Fei et al., 2023).(modified after HLNUG, 2007 andSeib et al., 2022).
During the late Cretaceous and Eocene, the Sprendlinger Horst was uplifted as a result of the Alpine Orogeny and the related rifting of the URG (Schumacher 2002;Derer et al., 2005;Dèzes et al., 2004).During the Miocene and Late Oligocene, the maximum horizontal stress shifted to an NNW orientation, thereby transtensive reactivation of NE-SW oriented fault systems occurred (Schumacher 2002).From the Miocene to recent times, the regional stress field is characterized by an NNW-SSE to NW-SE directed maximum horizontal stress leading to active subsidence and a transtensive regime in the northern URG (Schumacher 2002;Meixner et al., 2016).
Locally at the project site (Fig. 2b), the western part of the Lichtwiese Campus consists of granite and granodiorite of the Frankenstein Complex buried under a thin 2 -6 m cover of Quaternary sediments.The uppermost basement is characterized by a weathering zone with variable thicknesses of up to 30 m (Greifenhagen 2000).The eastern part of the Lichtwiese Campus consists of Permian sedimentary and volcanic rocks.The volcanic rocks are ophitic and vesicular andesites and basalts.The sedimentary units consist of alternating beds of arkoses, sandstones, and mudstones, which, however, are only found strongly weathered at the surface.Thus, no outcrops exist near the test site.To the northeast of

Methods
To minimize the shortcomings of individual exploration methods, emphasis is placed on combining independent techniques and results, to converge to a common interpretation of sub-surface structure and properties.Specifically, Electrical Resistivity Tomography (ERT) and gravimetric measurements were carried out to support a 2D seismic survey (Fig. 3).
Local outcrops were investigated as analogues to obtain petrophysical data and to develop the conceptual reservoir model.Additional lithostratigraphic information was obtained from a construction-related outcrop at the campus.

Geotechnical borehole data
400 archive records from local geotechnical drillings were manually analysed to infer the surface distribution of geological units (HLNUG, 2020;Seib et al., 2022).82% of the boreholes do not exceed a depth of 20 m, which only allows inferences about the near-surface geology.Locally, boreholes of up to 50 m depth are present (Fig. 4a).

Geothermal boreholes and surface outcrop
The drilling operations on the MD-BTES were performed from June to November 2022.They were accompanied by cutting sampling in 3 m intervals.Between well completion and installation of the inner pipe, borehole geophysical measurements and seismic tomography were carried out.In August 2023, the installation was completed with the installation of the inner pipe of the coaxial borehole heat exchanger.This study reports only the lithological results of the cutting samples.
By chance, an excavation pit of a university building construction was opened near the drill site shortly after the drilling.The pit allowed the sampling of material for lab measurements and geological characterization.

Petrophysical data
Petrophysical parameters are important for the interpretation of geophysical measurements and the parameterization of geological and numerical models.Additional data was extracted from extensive outcrop analogue studies in the northern Odenwald and Sprendlinger Horst (Bär et al., 2011;Bär et al., 2020;Weinert et al., 2021).In total, petrophysical data from 640 samples are available (Weinert et al., 2021) just for the Sprendlinger Horst Thirty samples of Permian basalts were collected from the temporary excavation outcrop on the campus.The sample material was then prepared and analysed for thermal conductivity, grain, and bulk density as well as porosity at the Technical University of Darmstadt laboratory following the workflow presented in detail in Weydt et al. (2021).In total, 30 cylindrical samples with a dimension of 40 mm in diameter and ~30 mm in length according to ASTM D4543-19 (2019) were at hand.Before the lab measurements, the samples were oven-dried at 105 • C for more than 24 h and subsequently stored in a desiccator at room temperature (21 • C).The thermal conductivity of the matrix was measured with the Thermal Conductivity Scanning (TCS) method (Popov et al., 1999).Each sample was measured 6 times (3 times on each of the two parallel end faces) to account for sample heterogeneity.Grain density was analysed using a helium gas pycnometer (AccuPyc 1330) according to ASTM D5550-14 (2014).Bulk density was measured using a powder pycnometer (GeoPyc 1360).Both grain and bulk density were measured five times per sample.The gas-effective porosity of the samples was then calculated.

Gravimetry
In February and March 2022, gravimetric data was acquired at 267 stations using a Scintrex CG-6 Autograv gravimeter (Fig. 5).Along the seismic profiles (Fig. 5) the gravity field was measured with a point spacing of 10 m.At the storage site, raster measurements were performed with a spacing of 20 m.In the other parts of the Lichtwiese Campus, the point spacing was about 25 to 150 m.Base measurements were taken three times per day to record the instrument drift.The L. Seib et al. measurements were carried out over 10 days so that a total of 30 measurements were taken at the base station.Some gravity measurements were affected by various noise sources, e.g.neighboring construction sites.Accordingly, all values with a standard deviation of more than 0.06 mGal were excluded from the evaluation.The position was determined using a TOPCON differential GPS device.The measurements were calibrated at a fixed point, resulting in horizontal and vertical accuracy generally in the millimeter range.In the forest areas south of the Lichtwiese Campus, the accuracy was significantly lower, therefore, the height of the respective point was compared and corrected with a highresolution digital elevation model.
A Bouguer correction of the raw data was carried out with a reference density of 2670 kg m -3 .The topographic correction was done using the software GSolve.The gravity effect of the building's basements was examined, which is below the measurement's standard deviation for all points.Together with 329 reference data points over Germany provided by the LIAG and the HVBG (Hessian State Office of Land Management and Geo-Information), a grid of the Bouguer anomaly was calculated with a nominal resolution of 10 m using a cubic spline interpolator.
Since the Bouguer anomalies at Lichtwiese Campus are dominated by long-wavelength signals from several kilometer-deep crustal regions, a high-pass filter with a cutoff wavelength of 1000 mwasused to enhance the short-wavelength component of the gravity field and to detect the boundaries of near-surface geologic units (Cordell 1979). .
The Bouguer anomalies were correlated with the seismic data by creating a gravity forward model using IGMAS+ software (Götze and Lahmeyer 1988;Anikiev et al., 2020).IGMAS+ is a numerical simulation tool developed for the interpretation and modeling of potential field data.The gravity field is modeled on two 2D profiles along the seismic lines.The seismic reflectors thereby serve as a basis for defining boundaries of geological units in the model.By interactively changing the densities and/or geometries of the geological bodies along the working planes, the calculated gravity field can be fitted to the observed gravity field data.

2D seismic survey
Seismic reflection imaging is a commonly used method to resolve questions for, e.g., engineering tasks (Steeples and Miller 1988), hydrocarbon exploration (Brown and Fisher 1980), or geothermal questions (Hartmann et al., 2015).Most studies apply compressional waves  (P-waves) for higher penetration depth and a superior signal-to-noise ratio for a substantial part of the reflection signal, i.e. outside the ground roll cone defined by the relatively slow surface waves (Burschil and Buness 2020).Generally, the method works excellently in sedimentary environments when impedance contrasts are significant, seismic signals are coherent, and surfaces are only moderately dipping.Seismic methods have also been used to explore crystalline environments successfully (Milkereit et al., 1996;Malehmir et al., 2021), but contrasts in such environments are challenging to characterize due to often complex and inhomogeneous structures and low impedance contrasts (Cheraghi et al., 2021).
In the framework of the investigation, two P-wave reflection seismic profiles were acquired in March 2021 (Fig. 5).The LIAG hydraulic vibrator MHV4P was used as a vertical source together with 384 receivers, connected with 16 Geometrics Geodes as seismographs.The receivers, single vertical geophones SM-6, were spaced in 2 m intervals on two crossing lines in a single spread.The vibrator excited two linear sweep signals from 20 Hz to 240 Hz over 12 s periods every 4 m plus additional offset source points outside the receiver spread.The Geodes recorded 14 s listening time by a 1 ms sampling rate.The geometry of offset source points and fixed receiver spread lead to two profiles of m length (Profile 1) and 450 m length (Profile 2).Regarding the target depth of up to 750 m, the profile lengths are short, causing severe edge effects in later processing steps.However, the accessibility of the area did not allow for longer profiles.Processing was performed in Landmark Seisspace/Promax, with the parameters summarised in Table 1.A Migration Velocity Analysis (MVA) workflow for processing comparable to the Pre-stack Depth Migration (PSDM) by T was used.This workflow shows a superior coherency of reflectors compared to dip moveout processing and post-stack migration.

Electrical resistivity tomography
ERT is a widely used geophysical survey technique for sensing lithology (e.g.clay) and pore-water content.ERT has been successfully used for the characterization of near-surface aquifer structures in the context of weathered magmatic rocks (Singhal and Gupta 2010;Rønning et al., 2014;Comte et al., 2012).In this study, ERT was applied to enhance the knowledge of fault and deformation zones, weathering thicknesses, and groundwater.
A total of 18 2D profiles were acquired (Fig. 5b).Four profiles were already measured in a previous study (Seib 2021).The measurements were carried out with an ARES II automatic multi-electrode resistivity system with a maximum of 92 electrodes with a spacing between 1 m and 3 m.The contact resistance was low, ranging from 300 to 800 Ω for most of the profiles, yielding a generally good signal-to-noise ratio.In the profile segment located in forests with porous soils and high organic content, contact resistances up to 3000 Ω were measured.The ERT profiles were aligned perpendicular to the proposed geologic boundaries, crossing the location of the future SKEWS wellbores.The best survey success is achieved by the dipole-dipole electrode configuration due to its high resolution and sensitivity to geologic detail, especially for subvertical features such as faults (Seaton and Burbey 2002;Dahlin and Zhou 2004).The apparent resistivity data was processed with pyGIMLi (Rücker et al., 2017).pyGIMLi was used to invert the apparent resistivity data from the field into two-dimensional resistivity subsurface models.

Evaluation of geotechnical borehole reports
In the western part of the Lichtwiese Campus, the results of the geotechnical boreholes consistently reveal granitic or granodioritic material overlain by a weathered zone with a thickness from 5 m to m.In the east, Permian sedimentary and volcanic rocks are described in the borehole reports.A lithologic change from granitoid to volcanic rocks in the vicinity of the SKEWS wells can be inferred from the drilling data (Fig. 4b).

SKEWS geothermal boreholes and excavation pits
The SKEWS boreholes (Bossennec et al., prep.)   (Fig. 6).The Permian unit is heterogeneous and consists of vesicular and high porosity basalts and andesites.Locally, thin beds with 1 m to 5 m thickness of Permian fine-grained sand and siltstone were encountered.
The upper part of the basement, from 200 m to 470 m, consists of reddish granite, which in some areas is strongly altered.This reddish granite has a similar mineral content and appearance to the granite that occurs to the east of Darmstadt near the town of Dieburg at the Mainzer Berg quarry and is therefore in the following named Dieburg granite (Stein et al., 2022).From 460 m to 750 m, gray granodioritic material is encountered, which is associated with the granodiorite in the west of Darmstadt.In both granite units, brownish coloured, clay-rich sections were encountered that indicate advanced alteration of the surrounding host rock.
In the excavation pit altered and strongly fractured granite in the west and fractured and altered basalts in the east were outcropping (Fig. 7).Both units are separated by a red, clay-rich contact zone, likely representing the fault core.The granitoid partly shows internal breccia structures, strongly suggesting the damage zone of a fault structure.

Petrophysical parameters
In Table 2 and Fig. 8, the results of the petrophysical measurements are presented in accordance with data compilations from former research (Bär et al., 2020;Weinert et al., 2021).The granitoids from the Sprendlinger Horst exhibit grain and bulk densities of 2725 kg m -3 and 2660 kg m -3, respectively, with matrix porosities of about 2.1% (Table 2).The average thermal conductivity is 2.71 W m − 1 K − 1 (Weinert et al., 2021).The basalts sampled feature a comparatively lower grain and bulk density of 2697 kg m -3 and 2410 kg m -3 , respectively, but show a higher porosity of 10.65%.The average thermal conductivity is 1.74 W m − 1 K − 1 (Fig. 8).
The basalt samples collected in the excavation pit (Fig. 7) are considerably weathered and altered.Hence, the sample quality was not feasible for direct analogue data determination.Additional reference values of massive, greyish microcrystalline Basalt (Donnersberg-formation) from the public P 3 database for petrophysical rock properties (Bär et al., 2020;Hesse 2011) were used to display average values for the SKEWS subsurface.The basalts of the Donnersberg Formation exhibit a grain and bulk density of 2745 kg m -3 and 2622 kg m -3 with a porosity of 3.54%.The average thermal conductivity is 1.58 W m − 1 K − 1 .

Gravimetric Bouguer anomalies
For the whole region mapped, the Bouguer anomaly ranges from 1 to 20 mGal (Fig. 9a).The gravity field is dominated by an NW-SE-oriented regional trend.These variations of the gravity field are defined by the regional structure of the URG.
The residual gravimetric fields calculated by the filter sharp anomalies.The high-pass filter field (Fig. 9b) shows positive anomalies E and W of the SKEWS project site, ranging from 0.2 mGal to 0.45 mGal (Fig. 9b).Between these two maxima is a gravity minimum of − 0.15 mGal that continues northward.In addition, there is a local gravity field maximum central in the Lichtwiese Campus with values exceeding 0.4 mGal.The surrounding area offers − 0.2 mGal (Fig. 9b).Anomalies that are not well constrained by measuring stations are caused by the cubic spline interpolation technique.At these locations, the lack of data points leads to significant uncertainties.

2D seismic profiles
The seismic profiles with detectable signals reach a depth of about 600 m.Below 300 to 350 m, the quality of the data does not allow an appropriate interpretation.Both profiles are influenced by edge effects, e.g.caused by migration, and should therefore only be viewed with caution.
In P1, the NE area is characterized by multiple reflectors, one at a depth of 40 m and two at a depth of 90 m and 140 m (Fig. 10b, blue marker).Each reflector terminates on unstructured seismic facies with lower amplitude above and below.Several faults are interpreted in the center of the profile, based on offsets in seismic reflectors.Towards profile meter 280, the above-described reflectors disappear.Instead, there is a zone characterized by short, sub-horizontal reflectors framed by two fault zones (Fig. 10b, red marker).The SW part of the seismic section (Fig. 10b, green marker) is characterized by a series of subhorizontal reflectors to a depth of 50 m.Below this, the reflectors are short and unstructured.In P2 the SE area is characterized by continuous, dipping reflectors with slight offsets (Fig. 10b, blue marker).Above these dipping structures is an area with low amplitudes and without concise reflectors with a thickness of 100 m that is thinning out towards the NW.The NW part of the section consists of comparatively noncontinuous and mostly horizontal reflectors with high amplitudes (Fig. 10b, green marker).Based on the surface mapping, the interpreted faults, and the structure of the seismic reflectors, the eastern area of the profiles is assigned to basalts, while granite is encountered in the west (Fig. 10c).

Electrical resistivity tomography
The ERT profiles provide feasible information about the electrical resistivity of the upper 40 m of the MD-BTES site.The resistivity of the profiles ranges from 30 to 400 Ωm.A choice of a large number of profiles measured is presented in this study (Fig. 11).
Profiles 1, 3, and 4 are located directly at the MD-BTES site (Fig. 11b,  c, d).They display a zonation of resistivities with higher values around 100 Ωm in the NE and SW and low values of 50 Ωm in the center.In profile 4, there is an additional increase in resistivity compared to profiles 1 and 3, to about 150 Ωm in the central area at a depth of 20 to m.
Profile 7 (Fig. 11e) is perpendicular to the profiles 1 to 4. It shows an increase in resistivity towards the southeast with values of about Ωm.In the central region, there are resistivities of 100 Ωm.In the northwest, there is a resistivity minimum of about 40 Ωm up to a depth of 20 m.
Profile 9 (Fig. 11f) is located in the south of the MD-BTES site.It has a different resistivity structure with stronger resistivity contrasts compared to profiles 1, 3, and 4. In the northwest, there is a resistivity anomaly of approximately 350 Ωm.The SE is characterized by a zone of resistivities between 150 and 200 Ωm from a depth of about 15 m.

Gravimetry
In the east of the Lichtwiese Campus (Fig. 9d), which is sufficiently covered by gravity stations, there is a Bouguer gravity high that can be correlated with the observed Permian basalts.This basalt has direct contact with the granite NW of the boreholes as observed in the excavation pit (Fig. 7).West of the SKEWS site, this granite is causing a gravity high.The gravity low in between represents a zone of intense weathering adjacent to the fault zone (Fig. 9d).It can be traced to the north, which is interpreted as the continuation of this linear element.To  the south, a continuation of the lineament cannot be observed which may also be because of a lack of data.Alternatively, the change in Bouguer anomaly could be associated with a transition to granite south of the project site.Central on the Lichtwiese Campus the gravity maximum is regarded as the visualization of a basaltic dyke.However, the thickness of the Quaternary sediments and especially of the granodioritic weathering zone varies from about 3 to 20 m according to geotechnical drilling evidence.This can be an alternative interpretation of the center anomaly.

Correlation of seismic and gravimetry
In this study, the lab density data from analog samples was used to correlate seismic profiles and the gravimetric mapping results.For correlation, the high-pass filtered Bouguer anomaly was plotted over the seismic profiles.In-situ matrix densities can vary significantly from laboratory data measured at ambient T/P-conditions which is caused by fracturing and alteration (Weydt et al., 2022).For the granitic units, a range of values between 2520 kg m -3 for the near-surface and 2640 kg m -3 for the deeper granitic unit was defined.The Permian basalt ranges from 2450 to 2610 kg m -3 .
The seismic results indicate a lithological boundary at meter 200 (P1, Fig. 12) and 140 (P2, Fig. 13).Drilling data and cutting analysis of the SKEWS wells also indicate that a fault zone was crossed.
P1 in the NE intersects a positive gravimetric anomaly.It is interpreted as a sequence of basalt flows, which are represented by continuous reflectors.The modeled density (Fig. 8) variations of the basalt are correlated with decreasing weathering rate with depth.The mediumdense granite (2520 -2640 kg m -3 ) outcropping in the excavation pit causes an increase in gravity anomaly (SW of P1).Even though a fault zone is presumed (P1), no significant gravity minimum in the center of the profile is detected.This is likely to be caused by an orientation effect, as P1 is subparallel to the fault and is thus largely located in the damage zone.This reduces the observed density contrast (Fig. 14).
P2 (Fig. 13) has a gravity minimum in the central part of the profile.This is attributed to the fault damage zone-related weathering.As well as in P1, basalt layers are modeled in the SE.The outcropping granite is represented as in P1.The low amplitude zone SE in P2 is triggered by weathering and alteration effects and leads to a lack of acoustic impedance contrasts.

Electrical resistivity tomography
From the ERT data, a weathering horizon of the encountered Permian basalt can be deduced (Fig. 15).The observed lateral changes in resistivity are primarily attributed to differences in weathering and pore-water content.Due to the consistently low resistivities over the whole depth range, fresh material is not expected to be encountered even at the maximum penetration depth of 40 m.This state of pervasive weathering has been confirmed in the SKEWS boreholes.The known surface information and the drilling results were used to interpret the ERT measurements.
Profile 9 (Fig. 11e) is located south of the project site.Several pockets of high resistivity are observed in the center of the profile, which are interpreted as areas of fresh rock beneath the clay-rich weathering layer.The maximum resistivities encountered are about 400 Ωm, slightly higher than in the areas next to the boreholes.Profiles 1 -4 and 9 (Fig. 11, Fig. 16) are assigned to different lithologies.The former is assumed to be located within the area of weathered basalt (Fig. 16, 4a,  4b) while the latter is assigned to the granitic unit due to its overall higher resistivity and geological mapping results (Fig. 16, 9a, 9b).

Geological interpretation
Understanding the division of the Sprendlinger Horst into tectonically separated fault blocks is key for the geological interpretation at the Lichtwiese campus (Klemm 1938;Fahlbusch 1975Fahlbusch , 1980;;Marell 1989;Kowalczyk 2001;Hoffmann 2006).In the northern URG, NE-SW striking Variscan structures and N-S directed tectonic elements existed in the Permian, likely representing reactivated Variscan shear zones (Müller 1997).With the E-W-directed extensional stress field from Stefan onward (Henk 1993), this caused normal faulting on the N-S striking tectonic elements.A final dissection of the area in fault blocks and reactivation of Variscan and Permian fault systems presumably occurred during the Tertiary, associated with the development of the URG.Main phases of crustal extension took place in the middle to upper Eocene (E-W) (Villemin et al., 1986;Dèzes et al., 2004) the Oligocene (WSW-ENE to WNW-ESE) and the Miocene (NW-SE) (Buchmann and Connolly 2007;Frey et al., 2022).In this predominantly extensional setting, pre-existing Variscan fault systems were likely reactivated.A strike-slip component is possible, but could not be evidenced.
Based on this premise, a local, NNW-SSE striking graben structure is suspected for the SKEWS site (Fig. 17, Fig. 18), in which the location of volcanic activity during the Permian is associated with Variscan weakness zones.This presumed graben explains the preservation of the basalt sheets and resulting high thickness of Permian deposits.Regional variations in the thickness of the Rotliegend deposits are known, but such local variations have not yet been described at the Sprendlinger Horst (Müller 1997).The drilling data reveals sedimentary intercalations in the basalts (Fig. 18), indicating several eruptive phases as described in the literature (Al-Malabeh and Kempe 2009).In the updated geological map (Fig. 17), the basalt occurrence at the Lichtwiese Campus is connected with the northern outcrop at the Botanical Garden Campus based on literature (Fahlbusch 1980).

Impact on the storage system
The thermophysical formation properties influence the storage efficiency of the MD-BTES.The explored fault zone can negatively influence the storage performance of the BTES as fracture intensity is often enhanced (Seib et al., 2022).Special interest lies in the damage zone, where an increased flow path density can be expected (Caine et al., 1996).Note, that depending on the local characteristics, faults and fault cores may also act as fluid barriers (Caine et al., 2010;Bense et al., 2013;Scibek et al., 2016).In such a case, the fault zone would primarily influence the storage through the variable petrography of the surrounding rocks, for example, reduced thermal conductivity due to a high clay  The petrophysical properties of the basalt with thermal conductivities of 1.58 W m -1 K -1 to 1.74 W m -1 K -1 are not ideal for geothermal heat storage.Reduced injection and extraction performance must be expected (Welsch et al., 2016).The geoelectric, seismic, and borehole data indicate intensive fracturing and alteration, which promotes fluid flow and adversely affects the storage reservoir.The granitic and granodioritic units underlying the basalts exhibit superior properties for heat storage, averaging 2.71 W m − 1 K − 1 .In the borehole and the excavation pit, intensive alteration and a dense fracture network are observed in that unit.Based on the available data in this study, no reliable statement can be made about groundwater flow at depths below 200 m.

Assessment of the applied methods
To minimize the financial risks of installing a MD-BTES and to ensure good performance, an integrated exploration approach, especially in heterogeneous crystalline basement, is important for planning the well design and storage dimensions.The combination of methods chosen here is not necessarily optimal for every geological scenario.Instead, the optimal strategy should be tailored for each location based on the local geology and the data available, for which the experience gained here can be a valuable resource.
The evaluation of available geological data followed by analysis of local outcrops marks the first step of an investigation.The obtained geological information is used to select suitable exploration methods based on the petrophysical rock properties together with other limiting factors like local topography.At the same time, petrophysical data is of great importance for model parametrisation and therefore simulation and design of the storage system.
In this study, the analysis of the existing geotechnical drilling data helped to assess the general near-surface lithology at the local site scale and determine the optimal location of the geophysical profiles.Information on depth below 10 m was limited which meant that the description of the subsurface was previously mostly based on nearsurface conditions.The data still provided considerable added value for supplementing information from geological maps.
When assessing the seismic data, several points need to be considered.The length of the geophone spread was limited because of permission and natural preservation issues.This led to boundary effects, which had to be taken into consideration during interpretation.The vertical and horizontal resolution is limited by the sweep used, the rockspecific wave velocity, and the depth (Brown 2011;Faleide et al., 2021).Statements about small-scale changes on the meter scale, as they can occur in crystalline and volcanic rocks, cannot be made in this context.Nevertheless, the seismic profiles were indispensable for imaging the subsurface structure.They allowed the distinction of different volcanic and magmatic units and the identification of fault zones.To improve seismic imaging, acoustic core logs should be integrated as soon as data is available (Elger et al., 2021).Using a 3D seismic survey would have been beneficial due to its superior signal-to-noise ratio, and avoidance of out-of-plane reflections.As a 3D seismic survey provides a 3D data volume, elements like faults and geologic boundaries can be traced much easier in contrast to 2D lines, which often pose difficulties with fault correlation, as data is only visualized on sectional planes.However, 3D seismic surveys are typically more expensive than 2D lines, so the return on investment for storage projects with the lowest possible budget must be weighed on a project-specific basis.
The ERT profiles proved to be well applicable for the characterization of the shallow subsurface, which was not well imaged by the seismic campaign.The results indicated the groundwater level and clayey sections for the upper aquifer which was confirmed during drilling work.For a better assessment of the results, sensitivity tests can be carried out in the future to further increase the reliability of the results.
The gravimetric survey was conducted because the encountered sedimentary, volcanic, and magmatic rocks were expected to have significantly different bulk densities.In combination with the seismic measurements, gravimetry allowed the regionalization of the lithologies identified in the seismic data and the extraction of regional trends and structures.There are limitations when conducting gravimetric campaigns, which was problematic when conducting gravimetric campaigns in the southern part of the survey area where the dense forest resulted in a lower vertical resolution of the GPS.
Since distinguishing between the crystalline rocks of the Frankenstein Complex and sedimentary and volcanic units of the Sprendlinger Horst was a major goal of the survey, differentiation could have been made with magnetic measurements based on the magnetic susceptibility contrast between igneous and sedimentary rocks (Gabriel et al., 2011).An efficient possibility would have been an aeromagnetic survey, a widely used geophysical method that allows more efficient mapping of the crustal magnetic field than ground-based methods (Isles and Rankin 2014;Frey et al., 2023).However, obtaining a flying permit can be difficult in the vicinity of built-up areas.Locally, potentially disruptive influences on magnetic measurements from the nearby campus and the adjacent railroad line have to be considered.A further possibility to enhance resolution in larger depth is the use of electromagnetic methods.Magnetotellurics for example is sensitive to differences in electrical conductivity.It is typically used to explore hydrothermal geothermal systems with highly conductive brines (Geiermann and Schill 2010).In the context of exploration for MD-BTES, it can be used to identify fault zones with increased clay content, and water-bearing fractured zones, and for differentiation between rock formations.For all magnetic methods, measurements can be influenced by magnetic sources of interference such as buildings or railroad lines which adds certain constraints.
The costs incurred are a decisive factor in the selection of the exploration program.Based on the experience of the SKEWS project, the minimal and most cost-effective exploration approach consists of an exploratory borehole in combination with 2D seismic.Thereby both near-surface as well as deeper structures of the subsurfaces are sufficiently imaged.To reduce ambiguity of the seismic interpretation additional geophysical methods can be integrated.

Conclusions
In this study, the subsurface in the vicinity of a novel MD-BTES was explored based on a multidisciplinary workflow including gravimetric and 2D seismic surveys, petrophysical lab measurements of outcrop analog samples as well as ERT profiles.The findings from the investigation should help to streamline the exploration work on subsequent systems.
A potential fault area in the vicinity of the storage system was already postulated in advance to the geophysical survey.In the course of this study, the geologic situation was characterized in considerably more detail.The gathered knowledge was invaluable for the optimization of the drilling method.Due to the expected increased hydraulic permeabilities, it was decided to change from a hydraulic down-the-hole hammer to the rotary method.This quick decision reduced delays and drilling cost increases during the drilling phase.
Several geological findings were made which are of great importance for local geothermal projects in the future and the simulation of the storage system.
• The 750 m deep geothermal boreholes indicate a sequence of basaltic lavas near the surface, underlain by granitic rocks.In addition, an excavation pit 80 m northwest of the boreholes was analysed, which confirmed the basalt-granite contact zone.• A gravity low in the area of the MD-BTES wells is interpreted to be induced by the fault zone.West and east of the MD-BTES, gravity highs are correlated with decreasing fracturing and alteration intensity.• The seismic profiles reveal a complex subsurface structure of granitic and basaltic units, which are separated from each other by an eastward dipping normal fault.• ERT measurements indicate low electrical resistivities, which are most likely the result of weathering in combination with groundwater infiltration and flow.The ERT measurements also emphasize the heterogeneity of the near-surface geology.• The petrolophysical measurements reveal that the near-surface sections, including the Permian basalts, feature less favourable rock properties for heat storage, with high porosities and low thermal conductivities.However, the granitic units show a high potential for heat storage with low porosities and higher thermal conductivities.
This research emphasizes the critical role of integrating near-surface geophysics and a multidisciplinary interpretation of datasets with geological data for planning constructing medium-deep borehole heat exchangers.Compared to shallow BTES systems, such an increase in Fig. 17.Local geological map of the study area modified from survey results and field mapping (modified after HLNUG, 2007 andSeib et al., 2022).Line A-B locates the geological cross-section (Fig. 18).drilling depth affects the risks.Thus, these installations demand a strategic exploration approach to minimize costs for this emerging technology.Various geophysical techniques adapted to the specific geological context are pivotal for the enhancement of the threedimensional understanding of the medium-deep storage site's structure and potential hydrodynamic behavior.In the case examined, the most advantageous strategy involves pairing an exploratory borehole with a 2D seismic survey.Should the seismic survey yield insufficient data, electromagnetic measurements can complement it to capture rock units structures at greater depths and by gravimetric and magnetic potential field methods for regionalisation.

Fig. 3 .
Fig. 3. Concept for the exploration of the SKEWS site.

Fig. 4 .
Fig. 4. a) Location of geotechnical drilling with the corresponding depth at Lichtwiese campus The outlines of the different basement lithologies were extracted from the geological map of Hesse (HLNUG, 2007).b) Overview map showing the drilling locations of previous geotechnical boreholes at the Lichtwiese Campus.The colors correspond to the geology encountered below the Quaternary cover (modified with additional data after Seib et al., 2022).

Fig. 5 .
Fig. 5. Position of the geophysical measurements: a) Gravity Stations, b) ERT, and location of seismic common midpoints.The outlines of the different basement lithologies were extracted from the geological map of Hesse (HLNUG, 2007).

Fig. 6 .
Fig.6.Simplified geologic profile of the three wells based on cutting samples taken in 3 m intervals and mud logging data.

Fig. 7 .
Fig. 7. Location of an excavation pit nearby exposing a fault core presumably intersecting the SKEWS site.The red points in the map are the locations of the SKEWS boreholes.

Fig. 8 .
Fig. 8. Boxplot of the petrophysical parameters.Data for the granitoids of the Sprendlinger Horst were obtained from Weinert et al. (2021), while the petrophysical data for the basaltic lavas of the Donnersberg Fm. were retrieved from Bär et al. (2020).

Fig. 9 .
Fig. 9. a) Regional trend of the Boguer anomaly, b) 1000 m high-pass filtered Boguer anomaly.The outlines of the different basement lithologies were extracted from the geological map of Hesse (HLNUG, 2007).

Fig. 10 .
Fig. 10.Seismic profile P1 and P2 with the position of the wells and the profile intersection line.The top of the figure is a seismic reference level of 180 m elevation.a) Migrated profile, b) Interpreted section with marked faults in black and seismic facies marked in green, red, and blue, c) geologic interpretation.

Fig. 11 .
Fig. 11.Selected 2D profiles of the performed ERT measurements after inversion with logarithmic color scaling (b to f).The overview map (a) shows the positioning of the profiles.

Fig. 12 .
Fig. 12. Correlation of seismic profile 1 with the gravity data.a) the comparison between the measured and the modeled gravimetric anomaly.b) is the seismic profile.

Fig. 13 .
Fig. 13.Correlation of seismic profile 2 with the acquired gravity data.a) shows the comparison between the measured and the modeled gravimetric anomaly.b) is the seismic profile.

Fig. 14 .
Fig. 14. 3D-view of profiles P1 and P2 with the main fault and secondary faults.

Fig. 15 .
Fig. 15.3D image of the ERT measurements near the SKEWS boreholes.

Fig. 16 .
Fig. 16.Interpretation of two representative ERT profiles.Profile 4 is located directly on the project site and captures the upper weathered layer of the Permian basalt.Profile 9 is interpreted to be located in the granitic units.

Fig. 18 .
Fig. 18.Conceptual, not to scale ENE-WSW section through project area based on the results of the survey.The orientation of the cross-section is displayed in Fig.17.
Fig. 18.Conceptual, not to scale ENE-WSW section through project area based on the results of the survey.The orientation of the cross-section is displayed in Fig.17.
Fig. 18.Conceptual, not to scale ENE-WSW section through project area based on the results of the survey.The orientation of the cross-section is displayed in Fig.17.

Table 1
Seismic survey processing steps and parameters.

Table 2
Petrophysical parameters from outcrop sample lab measurements compared with literature values.