Characterizing a decametre-scale granitic reservoir using GPR and seismic methods – A case study for preparing hydraulic stimulations

Ground-penetrating radar (GPR) and seismic imaging have proven to be important tools for the characterization of rock volumes. Both methods provide information about the physical rock mass properties and geology structures away from boreholes or tunnel walls. Here, we present the results from a geophysical characterization campaign that was conducted in preparation for a decametre-scale hydraulic 15 stimulation experiment in the crystalline rock volume at the Grimsel Test Site (Central Switzerland). For this characterization experiment, we used tunnel based GPR reflection imaging as well as seismic traveltime tomography to investigate the volumes between several tunnels and boreholes. The interpretation of the GPR data with respect to geological structures is based on the unmigrated and migrated images. For the tomographic analysis of the seismic first-arrival traveltime data, we inverted for an aniso20 tropic velocity model described by the Thomsen parameters v0, ε and δ to account for the rock mass foliation. Subsequently, the GPR and seismic images were interpreted in combination with the geological model of the test volume and the known in-situ stress states. We found that the ductile shear zones are clearly imaged by GPR and show an increase in seismic anisotropy due to a stronger foliation, while they are not visible in the P-wave ( v0) velocity model. Regions of decreased seismic p-wave velocity, 25 however, correlate with regions of high fracture density. For geophysical characterization of potential deep geothermal reservoirs, our results imply that wireline compatible borehole GPR should be considered for shear zone characterization, and that seismic anisotropy and velocity information are desirable to acquire in order to gain information about ductile shear zones and fracture density, respectively. https://doi.org/10.5194/se-2020-40 Preprint. Discussion started: 30 March 2020 c © Author(s) 2020. CC BY 4.0 License.

Ground-penetrating radar (GPR) and seismic imaging have proven to be important tools for the characterization of rock volumes. Both methods provide information about the physical rock mass properties and geology structures away from boreholes or tunnel walls. Here, we present the results from a geophysical characterization campaign that was conducted in preparation for a decametre-scale hydraulic 15 stimulation experiment in the crystalline rock volume at the Grimsel Test Site (Central Switzerland). For this characterization experiment, we used tunnel based GPR reflection imaging as well as seismic traveltime tomography to investigate the volumes between several tunnels and boreholes. The interpretation of the GPR data with respect to geological structures is based on the unmigrated and migrated images. For the tomographic analysis of the seismic first-arrival traveltime data, we inverted for an aniso-20 tropic velocity model described by the Thomsen parameters 0 , and to account for the rock mass foliation. Subsequently, the GPR and seismic images were interpreted in combination with the geological model of the test volume and the known in-situ stress states. We found that the ductile shear zones are clearly imaged by GPR and show an increase in seismic anisotropy due to a stronger foliation, while they are not visible in the P-wave ( 0 ) velocity model. Regions of decreased seismic p-wave velocity, 25 however, correlate with regions of high fracture density. For geophysical characterization of potential deep geothermal reservoirs, our results imply that wireline compatible borehole GPR should be considered for shear zone characterization, and that seismic anisotropy and velocity information are desirable to acquire in order to gain information about ductile shear zones and fracture density, respectively.

30
Crystalline rock has been identified as the key host rock for geothermal energy exploitations (Brown et al., 2012) . The detailed geophysical and geological characterization of crystalline rock volumes are important steps in planning and managing geothermal reservoirs (e.g., Schmelzbach et al., 2016) . Before exploiting geothermal energy, the hydraulic transmissivity within the reservoir needs to be enhanced to allow for sufficient fluid flow at depths where temperatures are adequately high. Within crystalline 35 rock, fractures are the most important fluid flow pathways. Thus, the goal is to enhance their transmissivities and connectivities by high pressure fluid injections (so-called hydraulic stimulations). To this end, the key property to be characterised in crystalline rock is the geometrical characteristics, i.e. location, orientation, density as well as fracture length and aperture. However, quantification of the latest two is challenging. The local fracture densities may also lead to local heterogeneities in the in-situ stress 40 field (e.g. Bell et al. 1992, Valley andEvans 2010), and consequently impact the rock mass response to hydraulic stimulations. Additionally, the rock mass anisotropy (i.e. of elasticity, strength, etc.) is of importance as it can influence the propagation direction of induced fractures .
Available geological information originates mostly from the mapping of rock outcrops and tunnel walls, and from borehole data (e.g., geophysical and geological core and borehole logging). Direct sampling is 45 thus restricted to accessible locations leading to highly fragmented data sets with limited spatial coverage. To improve the spatial coverage inter-and extrapolation between these highly fragmented data are required. While such inter-and extrapolation can be reliably performed in settings with high borehole density (e.g., Krietsch et al., 2018a) , it can lead to high uncertainties in situations with sparse control data and/or significant spatial heterogeneity (Wellmann et al., 2010, 2014, Wellmann andCaumon. 50 2018) .
Geophysical imaging can help to infer geological structures in the rock mass away from boreholes and tunnels. In particular, ground penetrating radar (GPR), as well as seismic imaging and tomography have been shown to work well for tracing geological structures in crystalline rock (Schmelzbach et al., 2007). The challenge of imaging structures within crystalline rock with geophysical methods is that the rock 55 mass appears typically quite homogeneous on large scales with respect to its mechanical properties (i.e., no layering or geological sequencing is observable). Nevertheless, fractures and shear zones cause local perturbations of these mechanical properties, and dominate the hydraulic flow field and particle transport. As fractures only change the rock properties very locally (e.g. due to hydraulic and mechanical fracture apertures in the sub-millimeter range), individual fractures may remain undetected by geo-60 physical methods that have a spatial resolution in the metre-range, and hence only allow resolving bulk properties representative for a certain volume.
GPR is one of the highest-resolution geophysical imaging method that may allow recording the reflections from individual fractures in crystalline rock. Applications of GPR methods offer the potential to constrain both the geometry and hydraulic properties of fractured rock formations. Pioneering tests of GPR reflection and ray-based transmission methods for characterizing shear zones, fractures and intrusions in granitic rock date back to the 1980s Olsson et al., 1987Olsson et al., , 1988Olsson et al., , 1992Sandberg et al., 1989) . Experiments conducted in underground laboratories at Äspö in Sweden and the Grimsel Test Site (GTS) in Switzerland demonstrated that shear zones and lamprophyre dykes could be clearly identified in GPR images  . In salt tracer tests at GTS, the tracer flow through 70 permeable fractures was visualized using GPR cross-hole amplitude tomography (Niva et al., 1988) . More recently, Dorn et al. (2012) have shown how water-bearing fractures in a granitic rock could be imaged using multi-offset single-hole and cross-hole GPR techniques. As demonstrated in laboratory experiments (Grégoire and Hollender, 2004) , even the fracture aperture can be estimated from GPR data if the fracture filling is known. 75 Besides GPR investigations, seismic methods offer excellent opportunities to characterize crystalline rock masses as found at GTS. In particular, transmission experiments, such as cross-hole tomography, can provide valuable information on the presence or absence of fracture zones (Maurer and Green, 1997;Vasco et al., 1996Vasco et al., , 1998 . Seismic investigations at the GTS also provided clear evidence that the seismic anisotropy of the intact rock cannot be neglected (Vasco et al., 1998), which complicates the 80 analysis of seismic data. Nevertheless, the anisotropy parameters determined in rock volumes of intact rock and rock mass can provide valuable diagnostics on fracture orientations (e.g., Boadu and Long, 1996) . Numerical experiments by Rubino et al. (2017) show the first evidence that also the fracture connectivity can significantly reduce the anisotropy of a rock mass.
Theoretical considerations combined with experimental data show that seismic velocity is influenced by 85 the fracture density (i.e., an increased fracture density leads to a decrease in seismic velocity), the length of fractures intersecting a ray path, the fraction of the fracture in contact and the fracture filling (Boadu and Long, 1996) . There are, however, other factors altering seismic velocity, such as mineralogical rock composition, pore space and its filling, and the in-situ stress field. The fidelity of seismic velocity as a proxy for fracture density is thus rock type and site specific. 90 The link between seismic velocity and in-situ stresses is evident from laboratory studies (e.g. Holt et al., 1996) . For field studies, seismic velocities have been used to estimate in-situ stresses before drilling into hydrocarbon reservoirs (Sayers et al. 2002). However, it is not generally possible to infer in-situ stresses from seismic velocities due to the abovementioned factors that influence velocity. In contrast, transient poro-elastic stress changes can be imaged using seismic time-lapse tomography, as all factors except 95 effective stresses remain the same (Doetsch et al., 2018b;Rivet et al., 2016;Schopper et al., 2020) .The study presented here was conducted within the framework of the In-situ Stimulation and Circulation (ISC) project Doetsch et al., 2018a) in the crystalline rocks at the Grimsel Test Site (GTS), in central Switzerland. The ISC project addresses key scientific questions related to hydraulic stimulation aspects of enhanced geothermal systems (EGS), including a detailed reservoir characteriza-100 tion before and after the stimulation treatments. The crystalline target rock mass was carefully characterized with respect to its geological, hydraulic, mechanical and geophysical properties. Here we present the results of the pre-stimulation geophysical characterization and interpret them in the context of the geology  and in-situ stress field (Krietsch et al., 2018b). Prior to the experiments, they supported the planning and design of the stimulation experiments, and later, help improving the 105 interpretation of the stimulation experiments.

Site description
The Grimsel Test Site (GTS) is hosted within the crystalline rocks of the Central Aar Massif in central Switzerland, and operated by the Swiss national cooperative for disposal of radioactive waste (Nagra). The GTS is an underground research facility with an overburden of ~480 m. The test volume of the ISC 110 project has a size of approximately 20 m x 20 m x 20 m and is located at the very southern end of the laboratory and accessible from three tunnels ( Figure 1). For this study measurements only along the tunnel walls of the AU-and VE-tunnel were conducted. Both tunnels were drilled with a tunnel boring machine (TBM), resulting in an excavation damage zone of 0.4 m to 1.5 m (Frieg et al., 2012) and smooth walls that easily allow measurements on the tunnel walls. In addition to the tunnels, two boreholes dedi-115 cated for the high-pressure fluid injections (INJ-holes) and four boreholes dedicated for geophysical monitoring (GEO-holes) were used for geophysical borehole measurements (see Table 1).  https://doi.org/10.5194/se-2020-40 Preprint. Discussion started: 30 March 2020 c Author(s) 2020. CC BY 4.0 License.

Geology
The ISC test volume is located at the lithological boundary between the Central Aar Granite (CAGr) 125 and the Grimsel Granodiorite (GrGr) (Keusen et al., 1989) . These lithologies are close to the mineralogical transition between granitic and granodioritic rocks, and similar in quartz content (Wenning et al., 2018) . The major mineralogical difference between CAGr and GrGr is their amount in sheet silicate minerals (biotite and white mica) (Wehrens et al., 2016) . The rocks of these lithologies intruded the crystalline crust in the post-Variscan (Schaltegger and Corfu, 1992) and were subjected to Alpine defor-130 mation. During the latest phase of differentiation of the plutonic bodies, NW-striking SE-dipping zones of weakness were initiated, along which aplitic dykes and lamprophyres intruded (Steck, 1968;Wehrens et al., 2016) . The rock mass was metamorphosed with peak-pressures of 6.5 kBar and temperatures around 450 °C during the Alpine Orogeny (Goncalves et al., 2012) .
Simultaneously to the metamorphosis, the rock mass was deformed firstly in a ductile manner and later 135 during exhumation in a brittle-ductile manner. Two shear zone orientations can be distinguished within the test volume. The older shear zone orientation (S1) strikes NE-SW with a dip towards SE formed in thrusting regime under ductile conditions (Wehrens et al., 2016) (Figure 2). Parallel to this S1-orientation a pervasive foliation characterized by a preferential alignment of sheet silicates formed in the rock mass. 140 The younger shear zone orientation (S3) strikes EW with a sub-vertical dip towards south. S3 shows adextral strike-slip movement. Wherever the orientation of the meta-basic dykes (i.e. metamorphosed lamprophyres) is aligned with the S3-direction, the shear zones localized within those dykes ( Figure 2). During the late phase of the S3-oriented shearing, brittle fractures formed, as well as milky quartz veins (Wehrens et al., 2016) . These fractures often show a biotite coating, which indicates that this brittle de-145 formation took place under green schist conditions when biotite is chemically stable. The uplift of the Aar massif in the late alpine stage induced some partly open tension joints (Steck, 1968;Wehrens et al., 2016) . They seem to be the youngest features within the experimental test volume. Globally, the strike of the majority of the geological structures range between EW and NE-SW ( Figure 2). Krietsch et al. (2019a) identified four S1 and two S3 shear zones that cross cut the ISC test volume. Ad-150 ditionally, a highly fractured zone (up to 20 fracture per borehole meter) was identified in the eastern part of the test volume between the two S3 shear zones ( Figure 1).

155
In addition to the geology, the in-situ stress state within the test volume was characterized prior to the stimulation experiments (Krietsch et al., 2018b) . Based on the stress measurements, locations of two different principal stress measurement tensors were estimated for the description of the stress field. One tensor was measured far away (more than 30 m) from the shear zones (referred to as the unperturbed stress field), and the other tensor was derived from measurements at about 5 m from the shear zones 160 (perturbed stress field). The unperturbed tensor has a maximum principal stress magnitude of ~14 MPa which plunges with 40° towards East. The perturbed stress tensor indicates a maximum principal stress magnitude of 13 MPa that points sub-horizontally towards Southeast (135°). The minimum principal stress component drops from 8.7 MPa far away from the shear zone to <3 MPa at the S3 shear zones ( Table 1). The detailed analysis of the in-situ stress measurements is described in Krietsch et al. 165 (2018b) . Additional information on the in-situ stress field and its heterogeneity were obtained during the 12 main hydraulic stimulation experiments of the ISC project Krietsch et al., 2020) .

Geophysical methods and data acquisition
With the aim of imaging geological structures within the rock volume of the ISC project, several geophysical characterization techniques were applied. GPR reflection imaging and seismic traveltime tomography were chosen due to the sensitivity of these methods to the location and orientation of shear zones, as well as to the mechanical rock properties, such as fracture density. Previous tests have shown 175 that the acquisition of high-quality data is achievable in the environment of the Grimsel Test Site (e.g., Niva et al., 1988, Vasco et al., 1988, Maurer and Green, 1997. Other techniques, such as electrical resistance tomography (ERT) or frequency-domain electromagnetic (EM) methods, were considered, but not found suitable. For both methods, the electrical and metallic installations at the GTS are a major obstacle, in particular the grounding cables through the tunnels, 180 which channel the injected or induced electrical current, limiting the ability to investigate the highly resistive rock volume. In addition, for the ERT, the electrode-rock-coupling is difficult, as the rock has a very high resistivity (>10.000 Ωm), which makes galvanic coupling and injection of a current very difficult. For GPR, the high resistivity of the rock is a major advantage, as it ensures minimum damping of the propagating signal. 185

Ground Penetrating Radar (GPR)
The aim of the GPR imaging was to improve the geometrical information of the geological model for the ISC test volume. Structures within the rock volume of interest are the S1 and S3 shear zones, as well as individual fractures. In this study, we concentrated on the shear zones, as they cut through the entire test volume and were the target of hydraulic stimulations within the Grimsel ISC project. In order to im-190 prove the knowledge of the positions of the shear zones in the geological model, we jointly interpreted (i) unmigrated images and forward modelled data (Section 4.1) and (ii) the migrated and time-to-depth converted images (Section 4.2). This joint interpretation approach allowed constraining the location of the shear zones in 3D even though the GPR data were acquired only along 2D profiles. The interpretation of the 2D migrated reflection images on its own may be affected by out-of-plane artifacts resulting 195 from the geometry of subsurface features.

Data acquisition
GPR reflection data were acquired along the tunnel walls using shielded Malå HDR antennas with a nominal frequency of 160 MHz. The shielding of the antennas is important in this tunnel environment to ensure that strong interfering reflections from surface installations within the tunnels are attenuated. 200 Other shielded antennas with frequencies between 80 MHz and 250 MHz were tested as well, but the Malå HDR 160 MHz antennas were found to provide the best trade-off between resolution and depth penetration. Common-offset measurements with an antenna spacing of 0.33 m were performed in the VE and AU tunnels ( Figure 3a) with the imaging plane oriented in the direction of the experimental volume (approximately 45 from the vertical towards East and West, respectively). Traces were recorded 205 every 5 cm along the two approximately 50-m-long profiles. Data quality is high, due to the low signal damping as a result of the high resistivity of the rock and due to relatively few reflectors in the target volume.
210 Figure 3: Acquisition geometry for the three data sets. GPR data were acquired using shielded antennas in the tunnels (a), the 2D seismic data were acquired between the AU and VE tunnels (b) and the 3D seismic data were acquired with both sources and receivers located in tunnels and boreholes (c).

Processing of GPR data
Due to the very low electrical conductivity of the granitic rock, attenuation of the GPR signal was very 215 low and data quality high. Processing of the constant-offset GPR data was performed using Seis-Space/Promax® including the following steps (e.g., Schmelzbach et al.,   The workflow and its parameters were fine-tuned to enhance reflections in the GPR images. An initial data-independent amplitude scaling (in contrast to data-dependent scaling such as inverse envelope scal-230 ing) helped increasing signal amplitudes at later times, while not overly boosting remnants of surface reflections. The application of spectral whitening significantly suppresses mono-frequent ringing, whereas f-x deconvolution enhances coherent reflection events at the expense of random noise. A GPR velocity of 120 m/µs was used for the migration and time-to-depth conversion was derived from crosshole tomography and confirmed by testing different migration velocities by inspecting collapsing dif-235 fraction hyperbolae and the quality of other features. For the interpretation, the migrated images are used to infer he position and characteristics of the shear zones. Besides, the un-migrated data (processed up to step 8) are compared with modelled GPR data using the 3D geological model. The measurements on tunnel walls at unconventional (45° side-looking) angles, as well as the orientation of the shear zones and other 3D structures within the rock volume (tun-240 nels, boreholes) render the interpretation of the acquired GPR data challenging. Therefore, to verify reflections from known structures and to potentially improve the geometrical information of these structures, we modelled the reflections of known structures based on the available 3D geological model. For each source-receiver combination and based on the existing 3D geological model, a straight-ray based modelling algorithm (Schmelzbach et al., 2007) was used to calculate the 3D coordinates of the 245 reflection points on each reflector. Then, assuming a constant velocity of 120 m/µs, the reflection traveltimes were computed and a Ricker wavelet was placed at the time of the reflection. Under the assumption of constant velocity and planar reflectors, this modelling provides arrival times of reflected waves that can be compared with the observed (unmigrated) reflection traveltimes to infer on the location and orientation of the planar reflectors. Displaying the modelled reflection times on top of the processed (up 250 to step 8), but unmigrated data, allowed us to identify and verify various reflections from key geological structures. Any observed mismatch between the modelled and measured reflection traveltimes can be used to improve the 3D geological model.

Tunneltunnel seismics
Seismic data were recorded between the AU and the VE tunnels (Figure 3b). A total number of 120 100-Hz one-component geophones were installed at the VE tunnel walls in 0.5 m spacing, covering the western side of the experimental volume (Figure 3b). Seismic signals were generated in the AU tunnel using a small hammer and a chisel. The 120 source points were separated by 0.5 m and covered the 260 eastern side of the experiment. The high quality of recorded seismic data allowed the picking of 9500 first arrival traveltimes from the 14.400 total source-receiver combinations. The dominant frequency of the first-arriving seismic energy is ~1.1 kHz. We estimated the picking uncertainty to be around 0.04 ms, corresponding to roughly 0.5% of the traveltimes. Figure 4 shows a receiver gather for a geophone in the middle of the layout, but north 265 of the two S3 shear zones (location marked as "receiver 45" in Figure 3b). The receiver gather shows a clear delay in the arrival times for all shots fired south of the two S3 shear zones (shot numbers 75-100). Another delay can be observed for the shot positions 40-42, which are located at the intersection of shear zone S3.2 and the AU tunnel. These arrival delays give the first hint for influences of the shear zones on the seismic velocities. 270

275
Seismic data for the 3D tomography of the test volume were recorded using the 26 piezo-electric receivers installed for the passive seismic monitoring and sources in tunnels and boreholes (Figure 3c, , Schopper et al., 2020. A sparker source was used in the six water filled boreholes (INJ and GEO boreholes, see Table 1 and Figure 1b), resulting in 448 source positions with a source spacing of 0.5 m along the boreholes. In addition, eight hammers placed in the AU and VE tunnels were 280 used as sources. An example receiver gather (location marked as "receiver 5" in Figure 3c) for all sparker sources in the INJ1 borehole is shown in Figure 5. A total number of 10050 first arrival traveltimes were manually picked with an estimated uncertainty of 0.02 ms. In addition to the first arrivals that can be reliably picked, arrivals with a linear moveout in the receiver gather stand out in the data ( Figure 5). These events are recordings of secondary seismic waves that are generated by tube waves interacting  285 with open fractures intersecting the source borehole. The sparker source generates a strong pressure pulse in the water-filled borehole. Some of the energy is transmitted through the rock, but another large part of the energy travels as a pressure wave within the water of the borehole. As this pressure wave in the water interacts with open fractures intersecting the borehole, a part of the energy is transmitted into the surrounding rock. The fractures act as a source for secondary p-waves (tube waves converted to p-290 waves) that are recorded with a delay caused by the tube-wave travelling from the sparker position along the borehole to the intersection of the borehole with the fracture. The secondary p-waves have a linear-moveout, because the distance between the scatter point and the receiver remains constant and the tube-wave traveltime is a linear function of the source depth in the borehole. The apparent velocity of the linear-moveout is 1450 m/s, corresponding to the p-wave velocity in water. These linear-moveout 295 features can be used to identify and characterize borehole-intersecting fractures, as the amplitude of these secondary waves depends on the fracture aperture and compliance (Hunziker et al., 2019) .

Anisotropic 2D seismic traveltime tomography
Due to the foliation of the crystalline host rock within the test volume, the seismic velocity is anisotropic, which needs to be accounted for any traveltime tomography. For the 2D anisotropic traveltime 305 inversion, we extended the existing inversion framework presented by Doetsch et al. (2010) that uses the eikonal forward solver of Podvin and Lecomte (1991) and Tryggvason and Bergman (2006). We use the parameterization of Thomsen (1986) for weak seismic anisotropy in a transversely isotropic medium, with the velocity vp at an angle to the symmetry axis normal of the anisotropy plane given by: ( ) = 0 * (1 + 2 2 + 4 ).
(1) 310 0 is the velocity in the direction normal to the anisotropy plane, and and correspond to the strength of the anisotropy. represents the relative difference between the slow velocity 0 and the fast velocity (90°) within the anisotropy plane: (90°) 0 ⁄ = 1 + . The parameter is related to the ellipticity of the anisotropy. The inversion framework was set up in a way that 2D models of the parameters 0 , and were esti-315 mated simultaneously. It is also possible to simultaneously invert for the direction of the anisotropic symmetry axis, if this axis lies within the 2D inversion plane. If the anisotropic symmetry axis points out of the inversion plane, it is not possible to invert for the direction of this symmetry axis, because 2D data are not sufficient to constrain the two anisotropy angles. Therefore, in case of a 2D inversion, the direction of the symmetry axis has to be prescribed, and it is not possible to invert for it. For the Grim-320 sel test site, we know the direction of the seismic anisotropy from the direction of foliation, determined in the geological characterization. This direction, which is oblique to the inversion plane, was thus prescribed in the inversion. For the calculation of the ray paths, we assume that the ray paths are not affected by anisotropy. This assumption is only valid for crosshole or cross-tunnel geometries and if the anisotropy is weak, which is 325 given for our experiments. The traveltimes are calculated from the ray paths by calculating the traveltime in each model cell using Eq. 1 and summing these times up along each ray path. The sensitivity or Jacobian matrix for the anisotropy parameters are also calculated using Eq. 1, as well as the ray path information (see supplementary material for details). In addition to the estimation of the velocity and anisotropy parameters, the inversion algorithm also al-330 lows inverting for source and receiver static time shifts. These time shifts for individual sources and receivers can be compensate to local effects around the source and receiver positions, thereby improving the quality and consistency of the velocity and anisotropy parameter fields.

3D seismic traveltime tomography
For the anisotropic traveltime tomography of the seismic 3D data set, we perform a two-step inversion. 335 First we invert for global values of 0 , and , while prescribing the anisotropic symmetry axis. In a https://doi.org/10.5194/se-2020-40 Preprint. second step, we fix the values of and to the result of the first inversion in order to estimate the 3D heterogeneity of v0. We thus assume direction of the anisotropic symmetry axis to be known and the Thomsen parameters and (Eq. 1) to be constant throughout the test volume. These assumptions are necessary, as the 3D ray coverage in our experiment is not sufficient to invert for 3D distributions of 0 , 340 and as well as the two angles defining the symmetry axis. As for the 2D tomography, the symmetry axis was assumed to be parallel to foliation. In the first step of the inversion, and are estimated by calculating the velocity for each source-receiver pair by dividing the source-receiver distance by the traveltime and plotting these apparent velocities as a function of the angle between the ray path and the anisotropic symmetry axis. The assumed direction of the symmetry axis can thereby be verified by ana-345 lyzing the apparent velocity as a function of the angle to the symmetry axis. and are then estimated from fitting the anisotropic velocity of Eq. 1 to the angular velocity distribution. We found that this inversion problem is well constrained as it is highly over-determined, with three parameters ( , and 0 ) and a large number of traveltime observations (>10.000 in this study).
In the second inversion step, the values of , , as well as the direction of the symmetry axis are then 350 kept constant and only the heterogeneity in 0 is estimated. This simplification allows for algorithms built for isotropic 3D traveltime inversion to be used with minimal adjustment. Here, we use the inversion algorithm of Doetsch et al. (2010) that uses the eikonal solver of Podvin and Lecomte (1991) and Tryggvason and Bergman (2006) . 355 Figure 6: GPR reflection data measured from the AU tunnel in the N-S plane, looking west. These data are processed up to step 8 (Section 3.1.2) but not migrated. The shear zone locations inferred from the geological mapping are marked in the top panel. The bottom panel shows as an overlay GRP 360 reflections, which were modelled based on known shear zone locations and their average orientation shows the pre-processed (up to step 8), but un-migrated data acquired in the AU tunnel in the N-S plane, looking west. The three S1 shear zones are clearly visible in the radargram (arrows in top panel of Figure 6: GPR reflection data measured from the AU tunnel in the N-S plane, looking west. These data are processed up to step 8 (Section 3.1.2) but not migrated. The shear zone locations inferred from 365 the geological mapping are marked in the top panel. The bottom panel shows as an overlay GRP reflections, which were modelled based on known shear zone locations and their average orientation) and match the modelled reflections (lower panel) very well. The geological model was already of high quality and only needed minor updates of the shear-zone details to match the GPR data. Reflections from the two S3 shear zones are visible neither in the recorded nor in the modelled data, because these 370

Comparison of un-migrated and modelled GPR data
shear zones are perpendicular to the tunnel and thus do not create reflections for the side-looking GPR.

Migrated GPR data
The fully processed data acquired in the VE tunnel were used to analyse the more complex shear zone 380 geometry in this area, with the S1 shear zones being intersected by the two S3 shear zones. Figure 7: Fully processed and migrated GPR data acquired from the VE tunnel, along with its geological interpretation. Top and bottom panel show the same GPR image. Additionally, in the bottom panel, the intersections with the geological structures inferred from the geological model are shown. Here, especially the complex intersection zone of the S1 and S3 shear zones is of interest. shows the migrated 385 and time-to-depth converted GPR data along with our interpretation. The S3 shear zones that intersect the volume and cut through the S1 shear zones are evident as areas of low reflectivity. The S3 shear zones intersect the VE tunnel at steep angles, which makes them poor reflectors in constant-offset GPR images recorded from the tunnel (compare with Fig. 6). Even the strongly fractured rock between the two S3 shear zones shows surprisingly little reflectivity, most likely due to fractures being parallel to 390 the S3 shear zones. Another reason for the little reflectivity might be the enrichment in phyllosilicates within the S3 shear zones, leading to an enhanced electrical conductivity (Wenning et al., 2018) . In contrast, borehole GPR data recorded from the GEO-boreholes (Figure 1 and Table 1) clearly image the S3 shear zones and even monitoring of saline tracer migration between them is possible using borehole GPR (Giertzuch et al., 2020) . 395 The S1 shear zones can be identified on both sides of the S3 shear zone in Figure 7, but cannot be traced through them. Also, S1-parallel brittle structures can be identified on both sides of the S3 shear zones. These types of S1-parallel structures of limited length (here 5-10 m) occur throughout the experiment volume and are also evident in borehole logs. These features are likely associated with brittle deformation along the foliation planes. During the time of retrograde brittle deformation of the S1 shear 400 zones, the S1 parallel foliation also might have been locally deformed under brittle conditions. The shear zone geometry as imaged by GPR fits with the previous geological model and adds detail in form of the S1-parallel features, which were not previously mapped. GPR also adds detail information in the area of the intersecting shear zones, which is not accessed by boreholes. 405

Anisotropic 2D seismic velocity model
The 9500 traveltimes (Section 3.2.1) were inverted for 2D distributions of the Thomsen parameters 0 , and (Section 3.2.3), while fixing the direction of the anisotropy axis to the direction of the foliation of the rock. Figure 8 shows the result of the inversion, with an RMS misfit of 0.04 ms (0.3% -0.6% of the traveltimes of 8-16 ms). In addition to the three parameter fields, source statics (standard deviation of 415 0.07 ms) were estimated to account for differences in the rock properties in the immediate vicinity of the source locations ( Figure 8d). Source statics are largest for sources within the shear zones (source numbers 40-42 and 70-75), where delays could already be observed in the raw waveforms (Figure 4), as well as in the AU cavern. In contrast to the tunnels, which were TBM drilled, the AU cavern was blasted, which leads to a larger excavation damage zone and thus larger source statics. 420 The p-wave velocity in direction of the anisotropic symmetry axis ( 0 ) in Figure 8a shows a low-velocity zone in East-West direction, adjacent to the AU tunnel as main feature. Otherwise, the velocity is relatively homogeneous (2%). The dominant low-velocity zone is located between the positions of the two S3 shear zones and extends about 15-20 m west from the AU tunnel. The NE to SW striking S1 shear zones are not visible in the 0 image. The model (Figure 8b) shows the relative increase of ve-425 locity within the anisotropy plane, compared to the velocity in direction of the symmetry axis ( 0 ). The model shows several high-anisotropy features in north-east to south-west direction. These features coincide in direction and location with the S1 shear zones. In addition to the previously known shear zone strands that connect the AU and VE tunnels, a high-feature can be observed just south of the location of the S3 shear zones near the AU tunnel. While undetected during the geological characterization 430  , this feature has recently been verified as a fourth S1 shear zone. Geologically, the S1 shear zones are characterized by an increase in foliation compared to the surrounding rock. As foliation is also the reason for the seismic anisotropy of the rock, it is no surprise that the increased foliation within the shear zones are detected as an increase of anisotropy. The model of  (Figure 8c), which characterizes the ellipticity of the anisotropy, shows only a minor variation and does 435 not contribute to our interpretation of the results.

3D velocity model
For the 3D seismic traveltime data set, we first estimate the anisotropy parameters that are then assumed to be constant throughout the volume. The red dots in Figure 9 show apparent velocity (source-receiver distance divided by traveltime) as a function of angle from the symmetry axis. This velocity distribution 445 confirms that the direction of the anisotropy axis is parallel to foliation, because the maximum velocity is at a 90° angle to the symmetry axis (i.e., within the anisotropy plane) and the distribution is symmetric about the 90° direction. A different symmetry axis would result in a shifted or skewed distribution. Fitting the anisotropy model of Eq. 1 to the apparent velocity distribution results in estimated value of =0.065 and = 0.038 and the velocity distribution is shown as a green line in Figure 9. The value 450 is consistent with the 2D result (Figure 8b), while the is slightly lower than the 2D result ( Figure 8c). is not very well resolved and thus not interpreted. The scatter of the apparent velocity values around the estimated distribution shows the inhomogeneity of the rock, which we explain in our inversion by variations in 0 . Keeping the and values fixed, we use the 10050 traveltimes to invert for the 3D distribution 0 . We use a 2 m x 2 m x 2 m grid size and fit the data to the estimated error level of 0.02 ms. 455 Figure 10a and b show cut sections through the 3D 0 model. Similar to the 2D image, the low-velocity zone between the S3 shear zones near the AU tunnel is the strongest feature. Horizontally, at the level of the tunnels, the low-velocity zone extends from the AU tunnel 15 m to the West, in agreement with the 2D results. Vertically, the low-velocity zone extends downwards near the AU tunnel and extends further to the West at a depth of 20 m below the tunnels (15 m elevation). While extending vertically 460 and in E-W direction, the N-S extent of the low-velocity zone is restricted to a few metersapproximately the extent between the two S3 shear zones. Similar to the 2D velocity model, there is no evidence of the S1 shear zones in the tomogram.

Integration of geo-mechanical a-priori knowledge and implications for hydraulic stimulations
The major aim of the GPR surveys was to constrain the knowledge on shear zones orientations and persistency for further spatial extrapolations away from tunnel and borehole intersections. The S1 shear zones could be traced in the unmigrated GPR image from the AU tunnel, as well as in the migrated image from the VE tunnel (see Sections 4.1 and 4.2). Additional S1 parallel fractures have been identified 480 inside the rock mass, which could not have been mapped during the geological characterization. This is in agreement with the observations from Krietsch et al. 2020 andVilliger et al. 2020. Both describe that fluid diffuses mostly along a single planar feature, during the hydraulic stimulation experiments targeting S1 structures. In contrast to this, the S3 shear zones are poorly visible in the GPR images. Above we named an increased fracture density between the S3 shear zones as potential explanation for the low re-485 flectivity. This would fit again the interpretations from Krietsch et al. (2020) and Villiger et al. (2019) , who argued that fluid diffuses within a complex fracture network, when injected into S3 shear zones. Thus, the GPR survey helped to refine the geological model published by Krietsch et al. (2018a) . From 2D seismic tomography, S1 and S3 could have been traced throughout the rock volume. The S1 shear zones are evident as an increase in anisotropy, visualized in the spatial variation of the Thomson 490 parameter . We argue that this is caused by the more distinct foliation within these shear zones, compared to the overall host rock. Due to this more distinct foliation, the rock anisotropy (i.e. the contrast in stiffness normal and parallel to the foliation) becomes larger. In contrast, the dense fracturing within the S3 shear zones reduces elasticity anisotropy and thus, are not visible in the variations. Nevertheless, they are clearly visible in the 0 tomogram as zones of reduced velocity. This velocity reduction can 495 have multiple reasons, for example changes in elastic properties of the rock mass and/or variation within the in-situ stress field. Wenning et al. (2018) described an enrichment in phyllosilicates as S3 shear zones are intersected. This, in combination with the enhanced fracture density between the S3 shear zones, may have resulted in locally reduced rock mass stiffness, strength and increased rock porosity, leading to a seismic velocity reduction. Interestingly, the zone of highest fracture density collo-500 cates with the lowest magnitudes. This is in agreement with Rubino et al. (2017) , who described that with increasing fracture connectivity the velocity anisotropy reduces.
To describe the influence of the rock mass elasticity, we compare the rheology and fracture density with the 0 measurements ( Figure 11). We extract the velocities along the boreholes from the 2D (borehole SBH4) and 3D (boreholes INJ1 and INJ2) tomograms and use the velocity 0 in the direction of the ani-505 sotropy axis, which coincides with the direction of foliation. Fractures were mapped along the boreholes using a combination of optical borehole televiewer images and core logs  . The fracture density values mapped along boreholes are biased by the borehole orientation with respect to fracture orientations, since e.g. borehole-parallel fractures are very likely to be missed in the mapping.
To overcome this influence, we transfer the longitudinal borehole observations (so-called P10) into volu-510 metric estimates (so-called P32) (Dershowitz and Herda, 1992) . Here, we follow the approach used by Brixel et al. (2020) for the same data set. Fracture frequency is generally low, mostly with three or fewer fractures per meter, but increases towards the shear zones. Along borehole SBH4 fracture density strongly increases as the S3 shear zones are penetrated. At a distance of 5 m from the shear zones along SBH4 0 starts to decrease, reaching its minimum between the 515 two shear zones. It is expected that seismic velocity decreases with increasing fracture density (e.g., Boadu and Long, 1996) . While our observations fit this expectation, also other factors influence seismic velocity. Thus, we cannot purely attribute the drop in 0 to the increase in fracture density. For example, the S3 shear zones that bind the zone of high fracture density consists of metabasic dykes, which include an increase in phyllosilicates, highly sheared material and fault gouge. As mentioned earlier, the 520 contrast in material compliance between host rock and metabasic dyke could have lalso ed to a reduction in 0 . Generally, seismic velocities in the tomograms are varying smoothly, which is primarily due to resolution limitations of the method as well as parameterization and regularization of the inversion. In contrast, measurements of fracture density depend on the aggregation interval and offset, but any two measurements are fully independent. Figures 10c and 11c show that all zones with low velocity are 525 characterized by increased fracture density. The main geological features for low seismic velocities are the S3 shear zones, which is characterized by a high fracture density in all boreholes. Additionally, the interval of 8-12 m in borehole INJ1 shows intermediate to low velocities and is characterized by above-average fracture density. At the same time, there are borehole locations that show a high fracture density, but the measured velocity is also high. 530 Nevertheless, since the highest fracture density in the rock mass was mapped where the 0 minimum was observed in the tomogram (Figure 11a), it seems valid to assume that the fracture density had a large impact on 0 at this location. At other locations between the S3 shear zones the fracture density did not increase so strongly. This fits the higher 0 observed between the shear zones in the tomogram (Figure 11a). However, 0 is reduced between and along the S3 shear zones compared to the host rock, 535 indicating the influence of the shear zones themselves. In contrast, we observe elevated fracture density at 43 m in INJ1 (Figure 10 c), that is not linked with a decrease in 0 . Here, we argue that the grid size for the 3D inversion of the seismic data is too coarse to capture local fracture density increases, which are not collocated with compliance contrasts. Our data thus supports the hypothesis that in our test volume, all zones with low seismic velocity have an above average fracture density, but not all areas of 540 high fracture density are images as a low-velocity zone. A direct interpretation of seismic velocities in terms of fracture density is thus not possible, but identified zones of low velocity are likely to have high fracture density. Villiger et al. (2019) stated that the stimulations targeting S3 shear zones are more seismogenic than injections into S1 shear zones and suggest that the geological setting of the S3 shear zones is decisive for it. We here argue too that this increased seismogenic behavior might be linked to the en-545 hanced fracture density and connectivity between the two S3 shear zones. Thus, we may also attribute an increased seismogenic behavior during high-pressure fluid injections to the zone were seismic wave velocities and anisotropy was reduced, compared to the overall host rock.
The extensive stress measurements of Krietsch et al. (2018b) enables us to also compare seismic veloc-550 ity to the in-situ stresses. The minimum principal stress σ3, measured in multiple locations along SBH4, are consistent with the 0 estimates. Both parameters start to decrease at a distance of 5 m to the S3 shear zones. We argue that both parameters show a similar behaviour, when approaching the shear zones, as both are strongly linked to the elastic properties of the shear zones (e.g. Boadu andLong, 1996, Krietsch el al., 2018b). Therefore, we conclude that the seismic velocity and the in-situ stresses 555 are linked indirectly via the elastic host rock properties. In-depth modelling of the geomechanics and stress state based on the geological and geophysical observations is the subject of further studies.

Conclusions
We combined GPR and seismic data to characterize a granitic rock volume in preparation for hydraulic stimulations. The GPR survey allowed imaging up to 25 m from the tunnels into the rock mass. We 565 could resolve the S1 shear zones that consist of an enhanced transversal isotropic elasticity model, compared to the host rock. Shear zones and areas with higher fracture density and connectivity were characterized by poor GPR reflectivity. The analysis of the seismic data required an anisotropic inversion of the traveltime tomography, due to the pervasive foliation within the rock mass. The seismic tomography indicates that the zones of highest anisotropy are the S1 shear zones, which are also clearly visible in the 570 GPR images. These structures are rather planar and contain few distinct fractures. In contrast, the S3 shear zones are visualized by reduced seismic anisotropy, which might be linked to locally enhanced fracture connectivity. Similarly, the S3 shear zones are characterized by reduced seismic wave velocities with respect to the host rock. A comparison between the in-situ stress state and seismic velocity showed that stress decreases if seismic velocity decreases. Possibly, more intense fracturing is associ-575 ated with a reduction in compliance leading to reduced stresses. The results of the geophysical characterization might be linked to the observations from the hydraulic stimulation experiments, as seismic responses and the propagation behaviour of the injection fluid seem to correlate with the observed seismic velocity and anisotropy around the injection location. Therefore, we argue that the geophysical characterization of the ISC test volume using a combination of GPR and seismic methods provided val-580 uable insight into the rock mass geology and was important for the understanding of the hydraulic stimulation experiments. For geophysical characterization prior to development of potential deep geothermal reservoirs, we suggest to consider borehole GPR, which is available for deep boreholes, for imaging shear zones. Seismic surveys should be designed to resolve velocity as well as anisotropy variations so that fracture zones as well as ductile shear zones can be located. 585