Characterization of the highly fractured zone at the Grimsel Test Site based on hydraulic tomography

In this study, we infer the structural and hydraulic properties of the highly fractured zone at the Grimsel Test Site in Switzerland using a stochastic inversion method. The fractured rock is modeled directly as a discrete fracture network (DFN) within an impermeable rock matrix. Cross-hole transient pressure signals recorded from constant-rate injection tests at different intervals provide the basis for the (herein presented) first field application of the inversion. The experimental setup is realized by a multi-packer system. The geological mapping of the structures intercepted by boreholes as well as data from previous studies that were undertaken as part of the In Situ Stimulation and Circulation (ISC) experiments facilitate the setup of the site-dependent conceptual and forward model. The inversion results show that two preferential flow paths between the two boreholes can be distinguished: one is dominated by fractures with large hydraulic apertures, whereas the other path consists mainly of fractures with a smaller aperture. The probability of fractures linking both flow paths increases the closer we get to the second injection borehole. These results are in accordance with the findings of other studies conducted at the site during the ISC measurement campaign and add new insights into the highly fractured zone at this prominent study site.

1 Introduction

Solid rocks, such as in crystalline and bedrock formations, typically have a compact matrix of low permeability. Water pathways are focused on mechanical discontinuities that separate individual rock blocks over multiple scales. Such fractures are commonly described as planar structures and form a network that is hard to resolve at field sites. This is due to the high diversity and complexity of natural fracture networks, the difficulty involved with identifying fracture connectivities, and thus the difficulty involved with interpreting the hydraulic regime of an entire formation based on local fracture detection. Accordingly, fractured-aquifer characterization represents a challenge, with a relatively high cost related to the application of specialized field investigation techniques and to gathering a sufficient data set for reliable hydraulic description. The general poor understanding of how groundwater flows in fractured field sites is in contrast to the relevance of fractured environments that host elementary freshwater reservoirs worldwide (Chandra et al., 2019; Wilske et al., 2020; Spencer et al., 2021). Moreover, adequate characterization of the properties of fractured field sites concerns many subsurface engineering applications, such as the planning and operation of enhanced geothermal systems (Vogler et al., 2017; Kittilä et al., 2020), the evaluation of potential sites for a nuclear waste repositories (Follin et al., 2014; Li et al., 2022), or the description of an excavation-induced damaged zone around tunnels and openings (Armand et al., 2014; de La Vaissière et al., 2015).

Depending on the chosen experimental setting and the available data, different interpretations of the hydraulic and structural properties of a fracture network are possible. A fractured site can be inspected locally by borehole data (e.g., core mapping and geophysical image logs such as optical or acoustic televiewer). The depth and orientation of structures intercepted by boreholes characterize fracture intensity and prevalent fracture orientations (Armand et al., 2014; Krietsch et al., 2018; Chandra et al., 2019; Tan et al., 2020; Yin and Chen, 2020; Pavičić et al., 2021); furthermore, by fitting probability distributions to the parameters, a statistical analysis can be conducted (Barthélémy et al., 2009; Massiot et al., 2017). Single-hole and cross-hole flow and tracer tests are employed to infer permeability and connectivity between different borehole intervals (Le Borgne et al., 2006; Follin et al., 2014; de La Vaissière et al., 2015; de La Bernardie et al., 2018; Jalali et al., 2018; Brixel et al., 2020b, a; Tan et al., 2020; Li et al., 2022), the velocity distribution (Kang et al., 2015), or transport properties (Kittilä et al., 2019; Lee et al., 2019).

Detailed insight into the properties of flow paths between adjacent boreholes can be gained by tomographic methods. The principle of all tomographic methods is perturbing the investigated system (e.g., by an injection of fluid, a tracer, a thermal anomaly, or an electric current) and recording the response at nearby receivers. In particular, geophysical tomographic methods are applied for the characterization of the rock properties, the identification of fractured (in particular highly fractured) zones, and the monitoring of flow pathways (Deparis et al., 2008; Dorn et al., 2012; Robinson et al., 2016; Doetsch et al., 2020). This is frequently done in combination with hydrogeological methods (Day-Lewis et al., 2003; Chen et al., 2006; Dorn et al., 2013; Voorn et al., 2015; Giertzuch et al., 2021b, a). A comprehensive portrayal of geophysical methods for the investigation of fractured field sites and the potential target applications is given in Day-Lewis et al. (2017).

In contrast to geophysical exploration techniques, hydraulic, pneumatic, or tracer tomography is based on a fluid or tracer injection at a source well. The response is recorded at different adjacent boreholes at different depth intervals. In most cases, the pressure signals or tracer arrival curves are evaluated by a continuous hydraulic conductivity distribution based on an equivalent porous media (EPM) concept (Yeh and Liu, 2000; Illman et al., 2008, 2009; Sharmeen et al., 2012; Zha et al., 2015, 2016; Zhao and Illman, 2017; Dong et al., 2019; Zhao et al., 2019; Kittilä et al., 2020; Tiedeman and Barrash, 2020; Poduri et al., 2021; Zhao et al., 2021; Jiang et al., 2022; Liu et al., 2022). Thus, detected high-conductivity zones correspond to the locations of fractures or faults. Further insights into the fracture properties and improved results can be gained by particle tracking simulations (Tiedeman and Barrash, 2020), binary priors representing either fracture or matrix (Poduri et al., 2021), or by generating synthetic models with similar features to the field site (Zha et al., 2015). Geostatistical methods apply a stochastic EPM, and different realizations of the subsurface are evaluated (Park et al., 2004; Blessent et al., 2011; Wang et al., 2017). Here, different facies represent different levels of fractured or intact rock, for which hydraulic conductivities are calibrated. In contrast to the EPM approach, the properties of the fracture network are inferred more directly by calibrating a connectivity pattern (Fischer et al., 2018b, a; Klepikova et al., 2020).

Our inversion approach differs from previous studies insofar as the fractured rock is represented explicitly as a discrete fracture network (DFN) and the hydraulic and structural parameters of the fractures are inferred directly. The great number of unknown parameters prevents the minimization of an objective function between simulated and observed data, resulting in a single deterministic DFN. Instead, a stochastic approach is applied to consider the nonuniqueness of the results. This is accomplished by generating several realizations of the fracture network that are equally likely to be evaluated as a fracture probability map. The validity of the approach has been demonstrated for synthetic test cases in two dimensions (2D) (Somogyvári et al., 2017; Ringel et al., 2019) and three dimensions (3D) (Ringel et al., 2021). In this study, the new inversion method is applied to field data for the first time. We use transient pressure signals from hydraulic tomography experiments conducted as part of the In Situ Stimulation and Circulation (ISC) experiments at the Grimsel Test Site (GTS) in Switzerland. Proper evaluation and validation of a new approach requires controlled tests, and the GTS and ISC experiments pose a well-explored site for experimental validation. The objective of this paper is to reveal the feasibility and capability of 3D DFN inversion using a small-scale example. This study provides an elementary link between the theoretical development of a new inversion algorithm based on synthetic test cases and field applications, although the small scale may not be representative of the much larger scale of groundwater reservoirs.

The paper is structured as follows: in Sect. 2, we describe the site and the hydraulic tomography experiments to be used for the inversion; the implementation of the inversion is elaborated upon in Sect. 3; we then review the forward modeling procedure (Sect. 3.1) and the general inversion framework (Sect. 3.2) developed in previous works with synthetic test cases; in Sect. 3.3 and 3.4, we explain the site-dependent inversion setting (i.e., the conceptual model and the prior parameter distributions that serve as basis for a stochastic inversion procedure) and discuss and justify the necessary constraints and assumptions; finally, the inversion results are interpreted and compared with findings from related ISC experiments in Sect. 4.

Figure 1 (a) General overview of the ISC experimental site showing the AU (Auflockerungszone, i.e., excavation effects) tunnel, the VE (ventilation test) tunnel, all boreholes, and the two types of shear zones (Krietsch et al., 2018). (b) The volume that is investigated in this study, i.e., the zone between the two brittle–ductile (S3) faults.

2 Experimental setting

2.1 Test site

The GTS is an underground rock laboratory located in the Aar Massif in the Swiss Alps. The ISC experiments, which serve as the basis for this study, utilized 15 boreholes of 20–50 m depth, including two injection boreholes (Inj1 and Inj2). The other boreholes are used for stress and strain measurement as well as seismic, pressure, and temperature monitoring during the hydraulic stimulation phases (Krietsch et al., 2018). A general overview of the site showing the persistent structures and the boreholes is given in Fig. 1a. A summary of the experiments conducted during the ISC measurement campaign and their results are given in Amann et al. (2018) and Doetsch et al. (2018).

The crystalline rock in the southern part of the GTS (ISC experiment volume) has been moderately fractured. Ductile (S1) and brittle–ductile (S3) shear zones can be distinguished from the investigated rock volume (Fig. 1a; Krietsch et al., 2018). The shear zones consist of a fault core, a damage zone, and unperturbed host rock (Wenning et al., 2018). A 4–6 m highly fractured zone with a fracture density (P10) of around 3 m⁻¹ is present between the fault cores of the two S3 shear zones and is displayed in Fig. 1b.

The fractures can be distinguished in wall damage zones adjacent to the S3 faults and linking damage zones, i.e., fractures connecting both fault cores (Brixel et al., 2020b). Testing campaigns on the connectivity between several intervals of the injection boreholes revealed that the best response occurs between the intervals 3 and 4 of both injection boreholes, which are located in the aforementioned highly fractured zone. Therefore, this is not only a highly fractured zone but also the most permeable region with conductive fractures (Jalali et al., 2018). For this reason, the characterization of the hydraulic and structural properties of this region (Fig. 1b) is the target of this study. The geological mapping of the structures intercepted by the boreholes and tunnels provides the basis for the setup of the conceptual model (Krietsch et al., 2018).

2.2 Hydraulic tomography data

The hydraulic tomography tests that are applied in this study are part of the characterization phase of the ISC experiment. We utilize transient pressure signals from constant-rate injection tests in the intervals 3 and 4 of the injection boreholes Inj1 and Inj2. The different intervals are isolated by a multi-packer system. The properties of the packer intervals and the parameters of the injection are summarized in Table 1.

Table 1 Parameters of the packer intervals and the hydraulic tomography experiments.

Between each injection experiment, pressure recovery was possible. The pressure response of the fluid is measured using piezoresistive pressure transducers. The resolution of the pressure response data is approximately 0.5 kPa. The minimum principal stress is of the order of 8 MPa. As the injected fluid pressure is much below the minimum principal stress, the coupling between hydraulic and mechanical effects can be neglected in the forward modeling of the experiment. The fluid pressure is measured with Δt=2 s. In general, we use similar hydraulic tomography experiments as those applied by Klepikova et al. (2020) except for a shorter injection time (Table 1), which was chosen for computational reasons.

Figure 2 Pressure response in the different intervals provoked by a constant-rate injection applied sequentially to the intervals Inj1–Int3 (a), Inj1–Int4 (b), Inj2–Int3 (c), and Inj2–Int4 (d), according to Table 1. The pressure measured in the respective injection interval is shown on the left vertical axes, and the pressure signals measured in the observation intervals are given on the right vertical axes.

The pressure signals are shown in Fig. 2 for each injection interval. Due to the stochastic inversion approach, the noisy pressure response data can be directly utilized for the inversion without the necessity to smooth or filter the signals.

3 Implementation of the inversion

3.1 Forward modeling

Fractures are modeled as 2D objects with constant properties normal to the fracture midplane in a 3D rock matrix that is assumed to be impermeable. The pressure diffusion in a single fracture is described by

$$\begin{array}{}\text{(1)}& a\mathit{\rho }S{\displaystyle \frac{\partial p}{\partial t}}-{\mathrm{\nabla }}_T\cdot \left(a\mathit{\rho }{\displaystyle \frac{k_f}{\mathit{\mu }}}{\mathrm{\nabla }}_{T}p\right)=aq,\end{array}$$

where a (m) is the hydraulic aperture, ρ (kg m⁻³) is the density of the fluid, S (Pa⁻¹) is the specific storage, k_f (m²) is the permeability, μ (Pa s) is the dynamic viscosity, and q (kg m⁻³ s⁻¹) is a source/sink term. The pressure p (Pa) consists of the static pressure and the piezometric pressure. The permeability is related to the aperture by

$$\begin{array}{}\text{(2)}& k_f={\displaystyle \frac{a^{\mathrm{2}}}{\mathrm{12}}},\end{array}$$

and the subscript T of the gradient (∇_T) denotes that it is evaluated in the fracture plane (Zimmerman and Bodvarsson, 1996; Berre et al., 2019). In this study, flow in the shear zones is modeled using the same approach as flow in the DFN, i.e., the shear zones are represented as 2D objects whereby the flow parameters are given by hydraulic aperture and specific storage (Eq. 1). The equations are solved numerically using the finite element method (FEM) with a conforming discretization at the intersections of different fractures. The generation of the geometry and the meshing of the fractures and shear zones are implemented using the open-source mesh generator “Gmsh” (Geuzaine and Remacle, 2009). The geometry of each structure is created separately by the built-in geometry module of Gmsh. The fractures and the shear zones are connected for a conforming discretization at the intersections of different structures by the Boolean operations implemented in Gmsh.

The implemented boundary conditions are shown in Fig. 3 along with the S3 faults and the fractures intercepted by the injection boreholes obtained from optical televiewer logs (Krietsch et al., 2018).

Figure 3 Overview of the volume considered in the forward model and the boundary conditions (BCs). The geometry of the S3 faults is simplified to planes, and the fractures intercepted by the injection intervals are illustrated as plane ellipses.

The boundary conditions are chosen considering the fact that only a small volume of the ISC experiment is investigated in this study. Therefore, the following boundary conditions are applied:

• The AU (Auflockerungszone, i.e., excavation effects) tunnel is represented by a pressure boundary condition – in this case, ambient pressure.

• The way to the VE (ventilation test) tunnel cannot be modeled explicitly. Thus, we apply a Robin boundary condition as a transfer boundary condition to consider the transition of the flow and the extension of the shear zones towards the VE tunnel (Watanabe et al., 2017).

• A no-flow boundary condition is applied normal to the planes of the fractures and shear zones.

3.2 Inversion algorithm

The parameters of the DFN θ are treated as unknowns characterized by probability density functions. Based on the Bayesian equation, the posterior density function p(θ|d) of the parameters given the measured data d is proportional to the likelihood function

$$\begin{array}{}\text{(3)}& \log L\left(\mathit{d}\mathrm{|}\mathit{\theta }\right)\propto -\sum _{i=\mathrm{1}}^{N_{\mathrm{data}}}{\displaystyle \frac{{\left(d_i-f(\mathit{\theta }{)}_i\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}}}\end{array}$$

and the prior distributions p(θ) (Gelman et al., 2013). The term f(θ) refers to the forward simulation of the hydraulic tomography experiment for the DFN realization defined by the parameters θ. The posterior density function is evaluated by sampling from it according to the Markov chain Monte Carlo (MCMC) method. This is an iterative procedure whereby new samples θ^′ are proposed by a proposal function and accepted (${\mathit{\theta }}_i={\mathit{\theta }}^{\prime }$) with probability

$$\begin{array}{}\text{(4)}& \mathit{\alpha }=\min \left(\mathrm{1},{\displaystyle \frac{p\left({\mathit{\theta }}^{\prime }\mathrm{|}\mathit{d}\right)q\left({\mathit{\theta }}_{i-\mathrm{1}}\mathrm{|}{\mathit{\theta }}^{\prime }\right)}{p\left({\mathit{\theta }}_{i-\mathrm{1}}\mathrm{|}\mathit{d}\right)q\left({\mathit{\theta }}^{\prime }\mathrm{|}{\mathit{\theta }}_{i-\mathrm{1}}\right)}}\mathrm{|}\mathbf{J}\mathrm{|}\right)\end{array}$$

or rejected (${\mathit{\theta }}_i={\mathit{\theta }}_{i-\mathrm{1}}$) (Brooks et al., 2011). The so-called reversible jump MCMC algorithm allows one to change the number of parameters (Green, 1995). In this study, the number of parameters is adjusted by deleting or inserting a fracture within the prior bounds. The determinant of the Jacobian matrix |J| has to be considered for transdimensional updates. It equals 1 for parameters sampled from the prior without linking its value to preexisting parameters (Sambridge et al., 2006). The parameters of the inversion problem are adjusted by proposing a new value from a normal distribution whereby the mean of the normal distribution is given by the current value.

The variance σ² in the likelihood function (Eq. 3) accounts for different sources of uncertainties, such as measurement errors, modeling errors, and errors of the conceptual model. Therefore, the value of the variance is estimated separately for each pressure signal. This is implemented as part of the inversion algorithm after the update of the parameters of the DFN. The measured data are assumed to consist of a mean and a normally distributed error $\mathit{d}=\stackrel{\mathrm{‾}}{\mathit{d}}+\mathcal{N}(\mathrm{0},{\mathit{\sigma }}^{\mathrm{2}})$. With this assumption, the variance can be estimated by sampling from an inverse gamma distribution

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}^{\mathrm{2}}\mathrm{|}\mathit{d},\mathit{\theta }\sim \mathcal{I}\mathcal{G}\left({\displaystyle \frac{N_{\mathrm{data}}}{\mathrm{2}}},{\displaystyle \frac{\sum _{i=\mathrm{1}}^{N_{\mathrm{data}}}{\left(d_i-f(\mathit{\theta }{)}_i\right)}^{\mathrm{2}}}{\mathrm{2}}}\right),\end{array}$$

as introduced by Gelman (2006) and implemented by authors including Haario et al. (2006) and Ringel et al. (2019). For this reason, the noisy measured data can be directly utilized for the inversion without filtering or smoothing the signals.

In practice, one iteration of the inversion algorithm operates as follows: assuming that the insertion of a fracture is chosen in the MCMC algorithm, the parameters (position; length; fracture set, i.e., orientation; and hydraulic aperture) of the fracture are generated from the prior functions. The chosen parameters are evaluated by simulating the hydraulic tomography experiment with the proposed parameter set θ (i.e., including the new fracture). The outcome of the simulation is compared to the measured pressure signals. If the error is smaller (the likelihood, Eq. 3, is higher) or similar to the previous step (without the fracture), the acceptance probability (Eq. 4) is high (Ringel et al., 2021). After accepting or rejecting the proposed parameters, the variance is updated according to Eq. (5).

3.3 Inversion constraints

The overall inversion procedure relies on several simplifications concerning parameters with less importance for our research target. For instance, the parameters specifying the properties of the shear zones have to be fixed. In general, our aim is an optimal balance between the accuracy of the generated results and the computational cost of the inversion procedure.

The underlying conceptual model comprises simplifications of the properties of single fractures that serve as inversion constraints. We assume plane ellipses as the fracture shape, and the length of the minor axis equals half of the length of the major axis (i.e., the length ratio is fixed). The assumption of reducing the fracture shape to a 2D plane is a common assumption and is justified by the derivation of the cubic law and the large ratio between the fracture extensions and the fracture aperture (Zimmerman and Bodvarsson, 1996). The assumption of the fracture shape as an ellipse is reasonable because the flow is dominated by the path between the intersections of different fractures; therefore, no sharp edges are considered for the simulation of the flow in the DFN. The hydraulic aperture is assumed to be constant over the fracture plane. Two fracture sets are defined with fixed orientations based on the orientations of the structures intercepted by the two injection boreholes. Thus, the fracture set is chosen by the inversion algorithm for the fractures between the boreholes; however, the orientation assigned to the fracture sets is a default. Figure 4 shows the orientation of the structures between the S3 shear zones intercepted by the two injection boreholes and the orientations defined for the two fracture sets. The appearance and distribution of the fractures dominate the flow. Accordingly, the surrounding rock matrix is assumed to be impermeable.

Figure 4 Orientations of the structures between the fault cores of the S3 shear zones in the injection boreholes observed from optical televiewer logs (Krietsch et al., 2018), shown in gray and light blue, as well as the calculated orientations for the fracture sets applied for the conceptual model.

The investigated volume is limited to the volume between the two S3 shear zones (Fig. 1). The shear zones consist of a fault core and a damage zone. The permeability increases with distance from the fault core, where the cores are almost impermeable (Wenning et al., 2018). As the properties of the shear zones are not the target of this study, the shape is simplified and the associated hydraulic parameters are fixed. The shape of the shear zones is simplified to a plane rectangle (i.e., a linear interpolation between the shear zones' traces at the injection boreholes). A constant hydraulic aperture of $a_{\mathrm{SZ}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m is assigned. This small value is chosen based on preliminary in situ tests and the knowledge that the cores of the shear zones are impermeable at their tunnel intersection. A higher permeability of the shear zone at specific locations can be covered by placing fractures in the respective area that also accounts for the spatial variability in the permeability of the shear zone. Moreover, the specific storage value is fixed at $S_{\mathrm{SZ}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ Pa⁻¹. This high value is prescribed considering the results from cross-borehole tests (Klepikova et al., 2020). Fractures of fracture set 1 are approximately parallel to the S3 faults. Hence, a position close to an S3 fault also accounts for spatial changes in the permeability and specific storage of the S3 faults.

Overall, the application of constraints and assumptions about the fracture shape limit an exact reproduction of the structural properties of the tested rock mass. However, those parameters that have a major influence on the flow in the DFN are adjusted by the inversion algorithm within prescribed bounds. These are, in particular, the position and the hydraulic aperture of fractures. In contrast, parameters with minor effects on the flow behavior are fixed (e.g., the exact fracture orientation or the length ratio).

3.4 Prior distributions

The parameters to be inferred are the number of fractures, the position of the fractures, the fracture lengths, the respective hydraulic aperture for each fracture, and the specific storage coefficient that applies to the whole DFN. The specific storage S (Eq. 1) is given by the compressibility of water in theory (Freeze and Cherry, 1979). However, some fractures are only partially open; thus, due to the roughness of the surface, the specific storage can be increased compared with the theoretical value (Jalali et al., 2018). Moreover, the hydraulic aperture is generally smaller than the actual aperture (Berre et al., 2019). The specific storage is assumed to be valid for the whole DFN because two variable hydraulic parameters for each fracture are not feasible for the inversion algorithm. Accordingly, five different update types are implemented to be applied sequentially: the transdimensional update changes the number of parameters by either inserting a new fracture or deleting a fracture; the other update types keep the number of parameters constant but adjust position, length, hydraulic aperture, or the specific storage. For the update of the position, length, and hydraulic aperture, one fracture is chosen randomly, and a new value is proposed by a random perturbation of the current value.

Uniform prior distributions are applied (i.e., a parameter is specified by a constant probability between minimum and maximum possible values that are given in Table 2).

Table 2 Uniform prior distributions defined by a minimum and maximum possible value.

All spatial values refer to the position of the midpoint of the ellipse.

The spatial priors are derived in general from the position of fractures intersecting the injection boreholes. The maximum value in the x direction corresponds to the distance to the AU tunnel to apply the boundary condition. The prior for the north direction is given such that the fractures are located between the cores of the S3 shear zones. The elevation of fractures is expected to have a minor influence on the flow between the two boreholes, and a broader possible range for the elevation would be less resolved. In the following, x refers to easting+667 400 m, y refers to northing+158 800 m, and z refers to height+1700 m (Fig. 1). The minimum value for the fracture length is given by the borehole diameter, and the maximum possible value corresponds to the distance between the shear zones. Fractures proposed during iterative inversion which intersect with the fault cores of the shear zones are reduced to the part of the fracture within the investigated volume (Figs. 1b, 3).

The prior range for the aperture is approximated from the results of single- and cross-borehole tests (Jalali et al., 2018; Brixel et al., 2020b, a). The minimum specific storage value is given by the compressibility of water (Freeze and Cherry, 1979), whereas the maximum value is based on cross-borehole injection tests (Klepikova et al., 2020). Both prior distributions for the hydraulic parameters cause preferential flow in the DFN rather than in the shear zones, due to a smaller specific storage and a larger hydraulic aperture of the fractures.

Figure 5 Comparison of the observed pressure response with the simulation of the hydraulic tomography experiment for the posterior DFN realizations for injection in the intervals Inj1–Int3 (a), Inj1–Int4 (b), Inj2–Int3 (c), and Inj2–Int4 (d), according to Table 1. For visual clarity, the observed pressure signals have been smoothed.

4 Results

4.1 Processing of the results

Overall, 27 000 DFN realizations are considered to be posterior DFN realizations because they minimize the error and fulfill the prior conditions. DFN realizations from the initial 500 iterations are discarded as so-called “burn-in” iterations due to a higher error. The computation of the inversion was executed by an Intel Core i9 workstation with 10 cores and 128 GB RAM and took about 1 week. The posterior realizations are approximately equally likely. They reflect the uncertainty of the inversion results in contrast to a unique solution that would be obtained by a deterministic approach. To reduce the autocorrelation of the results, we keep every 100th realization for further processing, which is called “thinning” (Brooks et al., 2011). Using the stochastic approach applied here, the fit between the measured and simulated pressure signals of the hydraulic tomography experiment is evaluated by the posterior and prediction uncertainty. The posterior limits are calculated based on the simulated pressure signals of the posterior DFN realizations which correspond to the uncertainty of the inversion method. The uncertainty related to predicting new observations is a measure of the overall error as well as of conceptual simplifications, as it also considers the estimated variance (Eq. 5). The DFN realizations are evaluated using a fracture probability map (FPM) over the investigated volume. For this, the inspected rock volume is divided into raster elements. Each element records whether the element is part of a fracture. By taking the element-wise mean over all of the posterior DFN realizations, the probability that each raster element is part of a fracture is derived. The evaluation of the FPM summarizes the estimated position and length of the fractures (i.e., those parameters with major influence on the flow). The hydraulic aperture is evaluated on the same raster elements. If a raster element is part of a DFN realization, the respective aperture is taken from the DFN. Thus, the mean hydraulic aperture can be evaluated for each element.

4.2 Evaluation of the data

In the first step, the measured and simulated pressure signals are compared to assess the quality of the posterior realizations. Figure 5 shows the median fit and the 95 % limits of the forward simulation of the posterior DFN realizations and the 95 % limits of the prediction uncertainty along with the observed data.

Figure 5 demonstrates that the general shape and trend of the measured signals are reproduced by the simulated pressure curves checking the median fit and the 95 % posterior limits. This is especially the case for the response in interval 4 of both boreholes. The weaker fit of some signals in interval 3 indicates effects not covered by the inversion approach or forward simulations, such as deviations from the assumed fracture shape or fracture orientations. For a given DFN realization, the actual measured pressure signals are predicted. Due to measurement noise and simplifications concerning the DFN model, the 95 % limits of the prediction uncertainty are wider than for the posterior uncertainty.

Figure 6 Evaluation of the results from the fracture probability map (a) and the mean hydraulic aperture (b) for different cross sections (z). The boundaries of the investigated volume are indicated by the cuboid in the lower left.

4.3 Evaluation of the DFN realizations

The FPM and the mean hydraulic aperture are shown for different cross sections (z) in Fig. 6. The fractures intercepted by the injection intervals and the shear zones are fixed; therefore, they appear with a probability of 100 %. Their orientation, as derived from the optical televiewer logs, is assigned to these fractures; thus, the orientation is in the same range as the orientation defined for the fracture set, but the exact values vary.

Overall, two different connections with different levels of permeability are present. A flow path dominated by fractures with a large hydraulic aperture exists between injection interval 4 of both boreholes (Inj1–Int4 and Inj2–Int4). The fractures providing this connection are visible with a high probability in the cross sections z=16 m and mainly z=17 m. In general, a good respective permeable connection between two intervals is possible with a large hydraulic aperture of the fractures, with long fractures, with a long intersection length between different fractures, or with a correlation of these factors. In contrast, a connection with fractures with smaller hydraulic apertures appears between injection interval 3 of both boreholes (Inj1–Int3 and Inj2–Int3) and Inj2–Int4. This flow path is present with an average probability of approximately 50 % primarily in the cross section z=15 m. Fractures linking both flow paths appear more likely the closer the location is to injection borehole 2 (i.e., further east). The described behavior is also reflected in the measured data. All responses provoked by the injection in interval 4 of both boreholes are more distinct than for the injection in interval 3. Although a maximum hydraulic aperture of 10⁻³ m is enabled by the prior distribution, only a few fractures with a small probability appear with an aperture close to the maximum possible value, as visible in Fig. 6, at a depth of z=17 m. The specific storage coefficient converges to a mean value of $S=\mathrm{7.4}\times {\mathrm{10}}^{-\mathrm{7}}$ Pa⁻¹. Only a few updates were possible that occurred mainly during the burn-in iterations. Therefore, this value is interpreted as the result of an optimization (i.e., as the averaged specific storage to be applied for the whole DFN). The estimated specific storage is greater than the theoretical value that functioned as the minimum value of the prior distribution of the specific storage (Table 2). This considers a delay in the response that is not related exclusively to the compressibility of water (Freeze and Cherry, 1979) but also to, for example, the surface roughness or fractures that are only partially open. Multiplied by the maximum possible aperture (Table 2), the inferred value is well within the storativity range calculated from cross-borehole injection tests (Klepikova et al., 2020). Several fractures of fracture set 1 appear close to the S3.1 shear zone, indicating either permeable fractures close to the shear zone or a higher permeability of the shear zone in this region than the assigned value. This demonstrates that the prescribed assumptions with respect to the hydraulic properties of the shear zone do not induce crucial conceptual constraints in the inversion, but a locally high permeability of a shear zone is indicated by a locally high fracture probability.

Although the volume east of injection borehole 2 towards the AU tunnel is part of the inversion (i.e., fractures can be inserted or moved in this volume), the resolution of the results is low because various DFN realizations (i.e., fracture positions) are possible. Only the volume between the two injection boreholes can be evaluated with a sufficient resolution.

4.4 Comparison with other studies

The inferred flow paths consist of fractures with a high or rather low permeability, which is in accordance with the results of Klepikova et al. (2020). We also compare our results with the structures intersected by other boreholes drilled after conducting the experiments evaluated in this study. While this inversion approach is capable of identifying fractures that are hydraulically relevant, geophysical methods (such as optical televiewer logs) report all structures intercepted by boreholes independent of their permeability. The boreholes PRP1 and FBS1 are partially located within the prior range defined for this study. The 23–25 m depth interval for PRP1 has been identified as the interval with the highest transmissivity by Kittilä et al. (2019) and Brixel et al. (2020a). In 95 % of the posterior DFN realizations, at least one fracture is present in this interval. Fractures that intersect with the interval between the S3 faults of the FBS1 borehole are present in about 45 % of the posterior DFN realizations. This supports the fact that crucial hydraulic features of the DFN can be identified by the presented inversion approach. Still, even if such successful local validation is possible, there are no other independent measurements available to confirm the validity of the inverted complete DFN structure and its probability. Geophysical measurements, such as seismic data (Doetsch et al., 2020) or ground-penetrating radar (Giertzuch et al., 2021a), were able to characterize the ISC volume on a decameter scale and identify the persistent structures and the highly fractured zone; however, they could not delineate or specify the properties of single flow paths.

5 Conclusions

In this study, we characterized the highly fractured zone at the GTS based on transient pressure signals from hydraulic tomography experiments using a new stochastic inversion method. A stochastic approach was applied to assess the uncertainty of the measured data and the nonuniqueness of the results. The fractured rock is represented directly as a DFN model in the forward simulations. Several posterior DFN realizations that are approximately equally likely are evaluated, and two preferential flow paths dominated by a large or small hydraulic aperture are successfully identified. The presented method relies on some investigations that must be applied prior to the inversion (such as the mapping of structures intercepted by boreholes) and benefits from single- and cross-hole permeability tests for the definition of the hydraulic aperture range. If it is possible to further narrow down the prior range of the hydraulic parameters, the specific storage can be inferred separately for each fracture, instead of computing only a mean value for the whole DFN. In general, improved results and more insights into the fractured rock can be gained using the same inversion method but with more pressure signals from additional intervals and boreholes.

Future research is necessary on the generally most suitable definition of prior and proposal distributions, which are elementary for robust inversion and for deriving meaningful results. The efficiency of the MCMC sampling can be improved significantly by implementing more elaborate prior or proposal distributions, for example, relying on soft information and site-specific expertise. A further option is utilizing continuous inversion results (such as continuous hydraulic conductivity distributions) or geophysical measurements for highlighting a priori regions with a higher probability of the insertion of fractures or to define zones that are likely connected by fractures to reduce the number of necessary inversion iterations (Dong et al., 2019).

The introduced inversion framework can be applied in a highly flexible way for the characterization of different fractured sites by adapting the site-dependent parameters to meet the conditions of the tomography experiment at each site. Moreover, different types and sources of measured data can be processed for the inversion (such as tracer or in situ stress data), provided that a forward model is available that allows for the flexible update of DFN parameters. The workflow for the setup of the inversion problem is similar. The basis is the properties of the fractures intercepted by the boreholes, i.e., their position and orientation, obtained from optical or acoustic televiewer logs or outcrops. This knowledge is utilized for the prior distributions on the spatial parameters and for the specification of fracture sets. The prior distributions on the hydraulic parameters are based on cross-hole flow tests in this study. This can also be done by the evaluation of the hydraulic tomography experiments as a continuous hydraulic conductivity and specific storage tomogram. As the definition of priors and constraints delineates the range of feasible DFN realizations, this step has to be done carefully. However, the presented Bayesian framework allows the combination of multiple and diverse hard and soft data, which often exist in addition to hydraulic test data that are used to guide the inversion. As demonstrated here, overly tight constraints may be avoided by uniform prior distributions with large value ranges at the expense of a higher computational cost for the inversion. In practice, the amount of information describing the fractured rock is determined mainly by the hydraulic tomography data (i.e., by the number of intervals and boreholes).

The present study paves the way towards the applicability of the discrete inversion approach on a larger scale. The main issue will be to balance the degree of field testing with the desired fracture resolution and the associated computational cost. One possible direction is explicitly implementing only large conductive fractures. The role of smaller fractures with a lower permeability could be represented by calibrating a background permeability within the discrete fracture matrix approach (Berre et al., 2019). Another appealing direction is the representation of scale-dependent fracture sets by their statistical properties following a hierarchical parameterization (Ma et al., 2020).