Comparison of Hydraulic and Tracer Tomography for Discrete Fracture Network Inversion

Abstract

Fractures serve as highly conductive preferential flow paths for fluids in rocks, which are difficult to exactly reconstruct in numerical models. Especially, in low-conductive rocks, fractures are often the only pathways for advection of solutes and heat. The presented study compares the results from hydraulic and tracer tomography applied to invert a theoretical discrete fracture network (DFN) that is based on data from synthetic cross-well testing. For hydraulic tomography, pressure pulses in various injection intervals are induced and the pressure responses in the monitoring intervals of a nearby observation well are recorded. For tracer tomography, a conservative tracer is injected in different well levels and the depth-dependent breakthrough of the tracer is monitored. A recently introduced transdimensional Bayesian inversion procedure is applied for both tomographical methods, which adjusts the fracture positions, orientations, and numbers based on given geometrical fracture statistics. The used Metropolis-Hastings-Green algorithm is refined by the simultaneous estimation of the measurement error’s variance, that is, the measurement noise. Based on the presented application to invert the two-dimensional cross-section between source and the receiver well, the hydraulic tomography reveals itself to be more suitable for reconstructing the original DFN. This is based on a probabilistic representation of the inverted results by means of fracture probabilities.

1. Introduction

Tomographic methods are specifically suited for characterizing the spatial heterogeneity of aquifers. They are based on the inversion of signals that are recorded among multiple sources and receivers, which are employed to infer two-dimensional (2D) or three-dimensional (3D) images of the spatial hydraulic parameter distribution. For identifying those structural characteristics that are relevant for groundwater flow, solute, and heat transport, invasive concepts are available that apply hydraulic [1,2,3,4,5,6] or pneumatic pressure [7,8,9,10], or that use tracers, such as dye [11,12,13], salt [14,15,16,17], or heat tracers [18,19,20,21]. In comparison to many alternative near-surface geophysical applications, these signals are rarely introduced from the ground surface. Instead, pressure and tracer tomograms are typically recorded between the boreholes by cross-well testing. For example, Cardiff et al. [22], Bohling et al. [23], Brauchler et al. [2], and Paradis et al. [24] applied short-term pumping or slug tests in different well intervals or levels, which are isolated by packers. Recordings from pressure sensors in adjacent observation wells obtained tomographic data sets for two- or three-dimensional inversion of aquifer heterogeneity. In contrast, tracer tomography, as presented by Jiménez et al. [12] for solutes, and by Somogyvári and Bayer [25] for heat, is based on the multi-level injection of tracers. Here, the tracer breakthrough curves that were recorded in adjacent observation wells serve as input for spatial inversion procedures.

Most available tomographic inversion methods are especially suitable for the reconstruction of smoothly changing hydraulic properties, such as those that are often found in sedimentary deposits [1,4,6,26,27,28]. In comparison, fractured aquifers commonly exhibit sharp contrasts between highly conductive fractures and a low-permeable matrix. Tomograms that were generated with classical smoothing algorithms, such as the successive linear estimator, thus provide a blurred visualization of hydraulic properties, and they roughly estimate the fracture locations or connectivity [5,29,30,31,32,33,34]. As alternatives, a hybrid procedure was presented by Wang et al. [35], and a travel-time based approach by Brauchler et al. [36] for exploring the trend and location of highly conductive channels. However, all of these inversion techniques similarly pixelate the spatial hydraulic properties and thus offer limited insight into the geometric properties of fractures.

Due to that, we chose to rely on a DFN approach, since this method is more realistic when emphasizing the features of a DFN and the connectivity of the fracture network. In contrast, equivalent continuum models can be appropriated to characterize rocks with high fracture intensity. In this case, the approximation of the rock as highly heterogeneous porous medium holds. However, the DFN case study that is presented in this paper deals with a small number of fractures. Another significant advantage of the DFN approach is that this method allows for an easier observance of a priori known statistical information, such as fracture inclination and fracture intensity.

Recently, concepts have been published that focus on a discrete representation of inverted fractures. For example, Fischer et al. [37] introduced a method that is based on the use of the Cellular Automata concept for deterministic inversion of linear structures and applied it to systems with fractured and karstic conduits [38,39]. Klepikova et al. [21] calibrated discrete fracture networks (DFNs) with given horizontal and vertical fracture orientations, mainly for investigating fracture connectivity, but not for the adjustment of fracture network geometry. The method by Somogyvári et al. [40] iteratively tunes geometric fracture properties in order to match the tomographic tracer measurements. Fracture probability maps were generated based on synthetic DFN examples, which reveal the good reliability of inverted results for simple networks, but a minor capability for locating single fractures in complex networks. This demonstrates not only the highly demanding inversion problem when adjusting individual fractures in a network, but it also highlights the growing data demand of accurate tomographic inversion, with increasing network complexity.

Improved tomograms can be achieved by increasing the source-receiver density, by optimal and full usage of information of measured data, instead of, for example, utilizing the trends and characteristics of response data only, and by combining different measurements that deliver complementary data. The combined use of salt tracer and electrical resistivity tomography [17], of salt tracer and ground penetrating radar [41], as well as of nanoscale zero-valent iron tracer and flowmeter testing [42] have been proposed for better delineation of fracture flow paths. However, combining different methods also raises efforts and the cost of field investigation campaigns. Moreover, the comparison or coupling of different data inversions escalates the requirements for the tomographic inversion procedure.

The presented work builds upon the DFN inversion procedure by Somogyvári et al. [40]. The objective is to further develop the tomographic reconstruction of discrete fractures by comparing the results from the tracer-based inversion to that based on hydraulic pressure signals. Our inversion procedure is assigned to the stochastic methods. Therefore, the probabilistic characteristics of the results will be the same for different inversion runs. Nevertheless, individual results will change, although we use the same input data. In contrast, deterministic approaches (e.g., Klepikova et al. [21]) will produce the same results for every inversion run while using the same input parameters. In the following, the underlying Bayesian inversion procedure is first described and improved. Subsequently, a synthetic case study is introduced for testing the different tomographic concepts.

2. Methodology

The following section starts by introducing our case study, its geometrical parameters, and the settings of the experiment. The reconstruction of the DFN is done by transdimensional inversion, which is explained in the following part. We present a method that deals with noisy observation data at the end of this section.

2.1. DFN Case Study

A synthetic 2D DFN model is used for testing the inversion procedure. It is set up based on the fracture geometry of the Tschingelmad outcrop in the central Alps [43]. The simplified cross-section that is shown in Figure 1 provides different preferential flow paths between the given three source and three receiver points, and it hence represents an ideal example to test the specific capabilities of different tomographic methods. The sources on the left side may be interpreted as packed vertical intervals of a well, where the tracer or pressure signals are introduced. The receiver points on the right represent measurement points in an imaginary observation well. By systematic characterization of the fractures that are adjacent to the source and receiver wells and with the aid of packer tests, one could isolate the conductive and connective fractures and also realize such tests in practice. For the subsequent tomographic inversion, it is assumed that the profile in Figure 1 can be reproduced by the definition of two fracture sets. The sets differ with respect to the fracture orientation and aperture, as given in

Table 1. List of geometric constraints, experimental and inversion parameter settings (FLD: fracture length distribution).

Category	Parameter	Fracture Set #1	Fracture Set #2
Geometric constraints	Fracture inclination (°)	−19.48	74.48
	Fracture aperture (mm)	1.5	1
	FLD-mean (m)	9.9
	FLD-variance m2	8.5
Experimental parameter settings	Injection pressure (Pa)	3 × 105
	Injection concentration (mg/l)	40
Inversion parameter settings	Discretization length (m)	1
	$p_{add}/p_{del}/p_{shift}$	0.4/0.4/0.2
	Number of rjMCMC iterations	100,000

. These properties of the fracture sets and the fracture length distribution (FLD) are geometric constraints of the inversion. Within these constraints, the inversion procedure will explore the suitable fracture combinations.

Theoretical experiments are simulated with different source signals to obtain receiver data that is used for the inversion to calibrate the DFN. Table 1 sums up the details of these experiments. Hydraulic tomography is simulated with a series of independent pumping tests, using a Heaviside injection pressure signal (transient flow condition). Tracer tomography is simulated in a similar hydraulic regime as the hydraulic tomography case. Water is injected in the source wells using the same pressure gradient between the source and receiver wells as for the hydraulic tomography. We wait until the flow stabilizes after the injection of water, i.e., until the velocity field is stable (steady-state flow condition), before any tracer injection is applied. The tracer is continuously injected, with a constant concentration, until the termination of the experiment.

A fast and robust forward model is implemented to simulate the pressure signals and tracer breakthroughs. Darcy flow, which relates the pressure gradient to the velocity, and the continuity equation are used to calculate the pressure diffusion in the fracture network. The fracture aperture indicates the hydraulic conductivity due to Cubic law. The intersection between two or more fractures is considered by mass conservation. These equations are numerically solved by a finite difference approach with first-order accuracy in time and second-order accuracy in space [44,45]. We use an implicit time integration method, which is unconditionally stable. The advection-dispersion equation is numerically solved when considering the steady state velocity field. The numerical solution of the equation produces the concentration field of the experiment in the DFN. We obtain the tracer breakthroughs at the monitoring intervals in the receiver well by monitoring the concentration at specific points. The tracer transport is assumed to be conservative, and no loss towards the rock matrix is simulated. The focus of the present study is set on a conservative tracer, but the experiment could also be simulated and inverted with a heat tracer. In the model, heat conduction and heat convection are described by a partial differential equation that has the same form as the advection-diffusion equation. Additionally, however, the heat conduction from the water to the rock matrix needs to be taken into account.

2.2. Transdimensional Inversion

The DFN reconstruction is accomplished by using a transdimensional Markov Chain Monte Carlo methodology: the Metropolis-Hastings-Green algorithm (MHG), as introduced by Green [46]. The method is also known as reversible-jump Markov chain Monte Carlo (rjMCMC), as it operates with reversible model updates (jumps) that change the problem dimensions. This algorithm was implemented for DFN inversion by Somogyvári et al. [40].

rjMCMC uses an iterative structure with two main phases per iteration; an update and an evaluation phase. The flowchart in Figure 2 describes the inversion procedure and the workflow of our inversion algorithm. In the update phase, the parameters of the DFN model are modified by one of the following updates: fracture movement, fracture insertion, and fracture deletion. The DFN model parameter set includes the number of fractures and for each fracture its position, the fracture inclination, and the fracture length. The parameter set of a DFN realization is denoted by $\theta$. The dimension of the parameter space also changes, since the number of fractures varies through adding or deleting a fracture. This motivates the description of transdimensional inversion. We estimate the number of fractures during the iterations of the inversion to consider the computational costs of the inversion problem. The number of fractures is in the range of $10$ to ${10}^2$ and the dimension of the parameter space is three times the number of fractures. These parameters have to be adjusted by the inversion algorithm to reproduce the measured pressure signals or tracer breakthroughs.

The kind of update is randomly selected, with a predefined probability $p_{add}/p_{del}/p_{shift}$, which may change during the inversion process, e.g., if the fracture intensity is too high.

Fracture addition is the so-called birth move of the rjMCMC algorithm. For this update, first the DFN model gets discretized, adhering to a pre-defined spacing distance. This characteristic property of the fractured rock represents the shortest possible distance between two parallel fractures. A set of insertion points is identified based on this distance, where new fractures could be added. The insertion of fractures is only possible starting or intersecting with a given fracture since isolated fractures do not affect the transport. Counting the number of possible insertion points, we receive the insertion probability $P_{ins}$. The orientation, length, and physical properties of a new fracture are probabilistically chosen, based on the given characteristics of one or more fracture sets. In this study, a fracture set represents fractures with similar physical properties and inclinations (Table 1). The choice between different fracture sets is taken into account by the fracture set probability, $P_{set}$. However, further constraints and properties may be defined, depending on the information that is available for a specific case. The fracture length is drawn from the fracture length distribution (FLD) that is assumed for the corresponding fracture set, which gives us the probability of the length of a new fracture, $P_{length}$.

The point of intersection along the fracture (fracture position) is randomly selected after the geometric properties of the fracture are determined. The new fracture is inserted in the field by assigning the fracture position to a random insertion point. We obtain the probability of the position $P_{pos}$ where the new fracture intersects the existing DFN according to the number of possible points along the new fracture. The different probabilities, which have a part in the addition probability, are explained in detail in [40]. The resulting probability of fracture addition can be expressed as the product of the individual probabilities of the above-mentioned steps:

$$q_{add}=P_{set}{P}_{length}{P}_{ins}{P}_{pos}$$

(1)

Fracture deletion is the so-called death move of the rjMCMC algorithm. A randomly selected fracture gets deleted from the DFN model when it is applied. This update is constrained, only allowing the deletion of a fracture if the source and receiver points remain connected. This is required to properly run the forward model. The probability of this update reads as the reciprocal of the number of deletable fractures, $n_{deletable fractures}$:

$$q_{del}=\frac{1}{n_{deletable fractures}}$$

(2)

$q_{add}$ and $q_{del}$ are the transition probabilities. We receive the probability of the transition from the current DFN realization to a proposed realization by computing them. Therefore, $q_{add}$ and $q_{del}$ are the conditional probabilities of updating the DFN given the current realization.

Fracture movement is a straightforward update, representing a random perturbation of fracture coordinates. It does not change the dimension of the model parameter space. A fracture is randomly selected in our implementation, and is then shifted to another randomly drawn insertion point while using the same discretization as the fracture addition. Fracture movement updates symmetrically, that is, the probability of the forward and the reverse update is equal. The proposal ratio does not need to be calculated, because the fracture length and the fracture set will not change, i.e., only the likelihood ratio matters. In contrast, the parameters of a new fracture are drawn from the FLD and $P_{set}$ if we first delete a random fracture and add a new one in the next iteration. This requires the calculation of the proposal ratio for adding and deleting. Therefore, fracture moving necessitates less computational costs than adding or deleting. Hence, it is the most efficient update. This type of update is also known as the Metropolis update, as defined by Metropolis et al. [47].

The DFN model is iteratively adjusted in a loop that is based on the recorded pressure signals or tracer breakthroughs by one of the described updates, as it is shown in the flowchart (Figure 2). A proposed DFN model (denoted by the model parameters ${\theta }^{\prime }$) is evaluated while using the Metropolis-Hastings-Green acceptance criterion, i.e., the update probability $\alpha$, comparing the current DFN realization (denoted by $\theta$) with the proposed DFN realization:

$$\alpha =\min \left[1,\frac{P\left({\theta }^{\prime }\right)}{P\left(\theta \right)}\frac{L\left(\xi |{\theta }^{\prime }\right)}{L\left(\xi |\theta \right)}\frac{q\left({\theta }^{\prime }\to \theta \right)}{q\left(\theta \to {\theta }^{\prime }\right)}\left|J\right|\right].$$

(3)

Equation (3) is calculated from the product of the prior ratio $P\left({\theta }^{\prime }\right)/P\left(\theta \right)$, the likelihood ratio $L\left(\xi |{\theta }^{\prime }\right)/L(\xi |\theta )$ the probability of the data given the model parameters, the proposal ratio $q\left({\theta }^{\prime }\to \theta \right)/q\left(\theta \to {\theta }^{\prime }\right)$, and the determinant of the so-called Jacobian matrix $\left|J\right|$. The prior probability, multiplied by the likelihood function, yields the posterior probability, according to Bayes’ theorem. The Jacobian takes the change in the dimension of the parameter space into account, as provoked by fracture addition or fracture deletion. In the presented DFN inversion, the prior ratio is always equal to one, because, aside from the non-informative priors, all prior information is included in the proposal ratio. The likelihood ratio is calculated by dividing the likelihood function of the updated model by the non-updated model. The likelihood of a model can be expressed with the following Gaussian distribution if we assume that the observations are independent and follow a normal distribution:

$$L\left(\xi |\theta \right)=\frac{1}{\sqrt{4\pi {\sigma }^2}}\exp \left(-\frac{{\left({\xi }_{}-f\left(\theta \right)\right)}^2}{2{\sigma }^2}\right).$$

(4)

The transdimensional updates of the presented DFN inversion only deal with the discretized model parameter space. New fractures can only be inserted to discrete insertion points, with discretized fracture length distribution. Deleting a fracture from a DFN is also a discrete operation. The Jacobian matrix of these discrete updates is always equal to one [48]. This makes this inversion approach very efficient, as the transdimensional updates can be evaluated at the same computational cost as updates that do not change the dimension of the parameter space. For the same reason, the fracture spacing discretizes the originally continuous fracture length distribution, since discretized values ensure an easy computation of the Jacobian in the acceptance probability $\alpha$, see Equation (3).

The MHG criteria ensure that the Markov chain is stationary, and it samples from the posterior probability distribution of the inverse problem. In our case, the MHG criteria simplify to:

$$\alpha =\min \left[1,\frac{L\left(\xi |{\theta }^{\prime }\right)}{L\left(\xi |\theta \right)}\frac{q\left({\theta }^{\prime }\to \theta \right)}{q\left(\theta \to {\theta }^{\prime }\right)}\right],$$

(5)

which is the same as the acceptance ratio of the Metropolis–Hastings algorithm.

The accepted realization is stored and the Markov chain continues in step. In the case of rejection, the previous model gets restored and it is used for the start of the next iteration. The final result of the inversion is the ensemble of model realizations along the whole Markov chain. However, in practice, the first part of the chain (burn-in phase) is discarded, as these realizations are influenced by the chosen initial model until the chain converges [49].

The likelihood governs the convergence rate of the inversion. In this case, the control parameter of the inversion is a fixed variance ${\sigma }^2$ of the normal distribution in Equation (4). Choosing a smaller variance value could enhance the convergence rate, with the risk of trapping the algorithm early in the local minimum. We inspect the acceptance probability $\alpha$ to study the effect of the variance, ${\sigma }^2$, on the convergence rate in more detail (Equation (3)). For this, we test the different settings of the variance. Note that, by changing the variance, only the likelihood ratio $L_{ratio}$ is affected, while the proposal and the prior ratio remain constant. First, we consider the case that the error between the observed and the current DFN simulated data decreases through a proposed realization. This yields a likelihood ratio that is greater than one, irrespective of whether a high or a low variance is chosen, and it induces a high acceptance probability.

Second, we suppose the choice of a high variance. The limit of $\frac{RM{S}_{current}^2-RM{S}_{proposed}^2}{2{\sigma }^2}$, is zero, as ${\mathsf{\sigma }}^2$ approaches infinite,

$$\underset{{\sigma }^2\to \infty }{\lim }\frac{RM{S}_{current}^2-RM{S}_{proposed}^2}{2{\sigma }^2}=0,$$

(6)

and therefore, the limit of $L_{ratio}$ is one $\left(\underset{{\sigma }^2\to \infty }{\lim }{L}_{ratio}=1\right)$. Hence, the use of a high variance causes a high acceptance rate

$$acc=\frac{N_{accepted}}{N_{iterations}}$$

(7)

independent of whether RMS becomes smaller or larger through the proposed DFN.

In contrast, the lower the variance is, the lower the mean RMS after the burn-in phase. However, a variance that is too low increases the possibility of the algorithm getting stuck in a local minimum as $L_{ratio}$ approaches infinity for a decreasing error and zero for an increasing RMS while using small variance. If a higher variance is chosen, the mean RMS increases as well, but more proposed realizations will be accepted and the mixing of different accepted realizations is better.

2.3. Estimation of the Noise Variance

In addition to the reconstruction of the DFN, we can estimate the variance of the measurement error, i.e., the measurement noise. While this value is the main control parameter of the convergence of the rjMCMC algorithm, in most implementations its value is empirically chosen. For pressure measurements, this value is relatively easy to properly choose, but for solute and especially for temperature measurements, it becomes difficult due to numerous environmental factors and the strong influence of imperfect experimental conditions (e.g., nonuniform physical properties of the borehole, sampling inconsistency, tracer degradation). Hence, we apply the following methodology for variance sampling: The measured data ${\xi }_{obs}$ (obtained by $n$ measurement points) is composed of the breakthrough curves (BTC) or the pressure signals that would be observed without noise $\tilde{\xi }$ and a measurement error $\epsilon$. Therefore, the recorded data can be expressed as ${\xi }_{obs}=\tilde{\xi }+\epsilon$. ${\epsilon }_i$ ($i=1,\dots ,n$) is independent and normally distributed, ${\epsilon }_i\sim N\left(0, {\sigma }_{noise}^2\right)$. The observation data is also normally distributed ${\xi }_{obs,i}\sim N\left({\tilde{\xi }}_i,{\sigma }_{noise}^2\right)$ with data mean $\tilde{\xi }$ and noise variance ${\sigma }_{noise}^2$.

Our aim is to determine $\tilde{\xi }$ and estimate the variance of the noise. Therefore, we need a probability density function (pdf) for ${\sigma }_{noise}^2$ that depends on the observed data and the model parameters to draw the variance from. The posterior pdf of ${\sigma }_{noise}^2$ can be calculated while using Bayes’ theorem with the model parameters $\theta$ (DFN and inversion parameters) and the observation data ${\xi }_{obs}$

$$p\left(\theta |{\xi }_{obs}\right)=\frac{p\left({\xi }_{obs}|\theta \right)p\left(\theta \right)}{p\left({\xi }_{obs}\right)}.$$

(8)

The denominator is assumed as constant and therefore can be neglected. For reasons of clarity, ${\mathsf{\theta }}_{-{\mathsf{\sigma }}_{noise}^2}$ denotes the model parameters without ${\mathsf{\sigma }}_{noise}^2$. Due to Bayes’ theorem and a constant $p\left({\xi }_{obs}\right),$ the posterior pdf of ${\sigma }_{noise}^2$ can be expressed as

$$p\left({\sigma }_{noise}^2|{\xi }_{obs}, {\theta }_{-{\sigma }_{noise}^2}\right) \sim p\left({\xi }_{obs}|{\sigma }_{noise}^2, {\theta }_{-{\sigma }_{noise}^2}\right) p\left({\sigma }_{noise}^2\right)$$

(9)

from which ${\sigma }_{noise}^2$ can be sampled [50,51].

To determine the likelihood function $p\left({\xi }_{obs}|{\sigma }_{noise}^2, {\theta }_{-{\sigma }_{noise}^2}\right)$, we assume that the breakthrough curves or the pressure signals of the current DFN realization (denoted by $f\left(\theta \right)$) calculated with the forward model are the data mean $\tilde{\xi }$. This assumption is valid after the burn-in phase, because the error between the data mean and the simulated data from the proposed DFN realization is small enough ${\tilde{\xi }}_i-f{\left(\theta \right)}_i\approx 0$. Therefore, the probability density function of the observation data given the model parameters, i.e., the likelihood function, can be calculated with Equation (4). The estimated variance is higher than the variance of the noise at the beginning of the inversion due to the higher error occurring during the first iterations, i.e., the burn-in phase. The estimated noise variance is inserted into the likelihood function to replace the fixed variance. This is valid because of the assumption ${\xi }_{obs,i} \sim N\left({\tilde{\xi }}_i,{\sigma }_{noise}^2\right),$ which is expressed by the likelihood function. Thus, the variance in the likelihood function is updated by estimating the noise variance in every iteration of the rjMCMC loop.

The use of a noninformative prior as pdf of ${\sigma }_{noise}^2$ ensures that the data mainly affect the posterior distribution. In this study, an inverse gamma (IG) and a uniform (U) prior are applied

$$p\left({\sigma }_{noise}^2\right)\sim IG\left(a,b\right),$$

(10)

$$p\left({\sigma }_{noise}^2\right)\sim U\left(0,c\right),$$

(11)

with the so-called hyperparameters $a$, $b$, and $c$ ($a$, $b$, $c$ $>$ $0$). A noninformative prior is obtained by using small values for $a$ and high values for $b$ and $c$. Bayes’ theorem leads to the following posterior pdfs that are dependent on the prior

$$p\left({\sigma }_{noise}^2|{\xi }_{obs},{\theta }_{-{\sigma }_{noise}^2}\right)\sim IG\left(a+\frac{n}{2},b+\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\xi }_{obs,i}-f{\left(\theta \right)}_i\right)}^2}{2}\right),$$

(12)

$$p\left({\sigma }_{noise}^2|{\xi }_{obs},{\theta }_{-{\sigma }_{noise}^2}\right)\sim IG\left(\frac{n}{2}-1, \frac{{{\displaystyle \sum }}_{i=1}^n{\left({\xi }_{obs,i}-f{\left(\theta \right)}_i\right)}^2}{2}\right).$$

(13)

The calculation of the posterior pdf is illustrated in Equation (12) while using the IG prior and in Equation (13) using the uniform prior. The IG prior has to be carefully handled, since the pdf is improper as $a$ and $1/b$ approach zero. A pdf is improper if the integral of the distribution is infinity, which violates the assumption that the cumulative probability distribution approaches one or any positive value for non-normalized distributions. In this case, the improper IG prior leads to improper posterior distribution. The results of the two different priors will be compared in the following because of this potential shortcoming [50,51].

After a new DFN realization is accepted or rejected (as explained in Section 2.2), the variance is estimated from the posterior pdf while using the Gibbs sampling method illustrated in Figure 3. The acceptance probability $\alpha$ of a Gibbs sampler is always one. Therefore, every proposed variance will be accepted [49]. The variance sampling is implemented in every iteration of the inversion loop after the evaluation of a proposed update of the DFN model parameters (see Figure 2). This is illustrated by adding the sampled variances to a histogram, as shown in step 4 of Figure 3. The described procedure is presented for an increasing number of iterations.

A similar procedure is employed by Fearnhead [52] and Punskaya et al. [53] to analyze a speech signal in order to obtain the real speech data without noise. Gallagher et al. [54] used the rjMCMC method to identify the abrupt changes in geochemical data series due to environmental or climatic variations and to fit a constant regression function to the recorded data curves. The regression function is defined piecewise between two changepoints, which model the abrupt changes in the data series. Adding, deleting, or moving a changepoint estimates the number and position of the changepoints. They estimate the variance of the noise while using an acceptance criterion to accept or reject a new variance proposal. In Sambridge [55], an rjMCMC algorithm is applied to reconstruct sea level variations assuming that the time series and the measured data contain stochastic noise.

3. Results

In the following, we first rely on a fixed variance in the likelihood function to compare the results of hydraulic and tracer tomography. In the second part, the results of the estimation of the noise variance are presented.

3.1. Results of the DFN Inversion

The result of the DFN inversion is not a single calibrated DFN model. In every iteration of the MCMC loop, a DFN realization is calculated and stored. The set of DFN realizations is called the ensemble. This result contains information regarding the probabilistic properties of the inverse problem, such as uncertainty and bi- or multimodality. The probabilistic characteristics stay the same for different inversion runs, as our inversion method is a stochastic method. To visualize the ensemble, we convert it to a fracture probability map [40]. This is done via the rasterization of the DFN realizations and by taking the mean of these raster models over the ensemble. The fracture probability of one point, i.e., the probability of this position to contain a fracture, is the number of DFN realizations that identify a fracture at this point divided by the whole number of realizations. This representation is equivalent to taking the mean of the ensemble, which is a common method of visualization (see [12,56]).

Figure 4 shows the tomographic reconstruction of the hydraulic tomography experiment. The fracture locations are well resolved, and the locations of the fracture intersections with the higher fracture probability roughly match those of the synthetic model. Fractures at the sides of the model (connecting the sources and receivers with the DFN) are fixed, thus they appear with 100% probability. This assumption is acceptable since it is possible to detect the fracture location and orientation in boreholes based on core and image logs (e.g., optical or acoustical televiewer) (see [45,57]).

The fracture probability map only shows the fractures with a probability higher than $10\%$ to present the results clearly. In Figure 4, the general characteristics of the original DFN can be identified. We can locate the area where the main horizontal connection between the source and receiver takes place. In addition, the non-fractured positions can be detected. Still, some errors are visible in the reconstruction. The quality of the recognition is better there so that the fractures on the right half of the model appear with higher probability because observation points are only placed on the right side of the model. The position of the fractures, which connect the fractures of the left side of the model horizontally with the middle of the investigated area, is uniformly distributed between $y\approx 8m$ and $y\approx 22m$. Therefore, we can conclude that this area contains fractures, but we cannot determine the exact position of a single fracture. In contrast, the position of the fractures, which horizontally connect the middle of the area with the right side of the model, has a higher probability and it is better resolved. The variance of a fracture position is much smaller on the right half of the investigated area. The position of the fracture with its center point at $x\approx 30m$ and $y\approx 12m$ can be fit to a normal distribution, since the fractures at the left and right side of it appear symmetrically with a decreasing probability. This corresponds to a Gaussian pdf of the position of this fracture. The aligned dots occur due to the rasterization in the fracture probability, and their size depends on the accuracy of converting a continuous DFN realization to the grid points of the probability raster. They represent fractures with a very low probability.

Figure 5 shows the tomographic reconstruction of the tracer tomography experiment. We check the acceptance rate of the rjMCMC loop after the burn-in phase to ensure that both of the results are comparable, see Equation (7). The acceptance rate in both cases is about $acc\approx 4.2\%$. From this, we can conclude that the different results of the different tomographic concepts do not depend on the parameter settings of the inversion, and especially not on the variance in the likelihood function. The most obvious difference to the hydraulic tomography example is the larger variability of the results. In contrast to the results of the hydraulic tomography (Figure 4), single fractures, except for fractures from the injection and receiver points, are less visible, since no fractures occur with a high probability. Additionally, neither the left nor the right side is resolved well with a high probability. Still, from the result in Figure 5, we get an idea regarding the fractured and non-fractured areas and the main characteristics of the DFN. The gaps appearing between some of the pixels are a matter of rasterization similar to the aligned dots and they do not occur in the individual DFN realizations. This results from a raster grid that does not accord well with the DFN geometry.

3.2. Results Estimating the Noise Variance

In contrast to the previous section, we assume that the observed tracer breakthroughs and pressure signals contain stochastic noise. Therefore, we want not only to do the inversion of the DFN, but to also estimate the noise variance and the data mean, i.e., the theoretically measured data without noise. This is accomplished according to Section 2.3. A synthetic example for noisy tracer breakthroughs with ${\sigma }_{noise}=3$ is shown in Figure 6.

The data mean is iteratively estimated by the inversion of the DFN by applying the forward model to calculate the tracer breakthroughs of the current DFN realization. The variance is sampled from the IG posterior pdf in every iteration of the inversion loop. These samples are illustrated in a histogram that compares this with the posterior pdf of ${\sigma }_{noise}^2$ $p\left({\sigma }_{noise}^2|{\xi }_{obs}, {\theta }_{-{\sigma }_{noise}^2}\right)$ of the last iteration of the rjMCMC loop in Figure 7. The pdf and the histogram in Figure 6 are as well not normalized, i.e., a cumulative probability function would approach any positive value for an increasing variance due to the non-normalized posterior pdf (see Equation (9)) from which the variance is sampled. The samples converge to a stationary distribution, since the histogram of the variance samples fulfills the posterior pdf of the last step. Both variance histograms are similar, as we compare the two different priors. Hence, we can conclude that both priors are in fact noninformative, because the posterior pdf is mainly affected by the data and not by the choice of the prior pdf. In addition, the prior IG pdf does not lead to an improper subsequent pdf in this case.

The same procedure could be applied to the noisy pressure signals (Figure 8). The histogram of the variance samples and the posterior pdf of ${\sigma }_{noise}^2$ are shown in Figure 9. Additionally, in this case, both of the priors provide similar histograms and posterior pdfs of ${\sigma }_{noise}^2$.

4. Discussion and Conclusions

The presented study deals with an exemplary hypothetical case, which makes it difficult to generalize the results to other experimental setups or general conclusions. Still, the application of the two tomographic concepts allows for the comparison of characteristic differences. Hydraulic tomography provides better results than tracer tomography, since the variability of the result is smaller and some fractures occur with a higher probability, i.e., the obtained fracture probability map of the hydraulic tomography is better resolved. In this case, the pressure signals are better for the inversion, because the complexity of the simulation is smaller. For practical use, many extra factors would probably make the tracer tomography better, e.g., the measurement noise is smaller, and tracer breakthrough curves can be more easily detected. Both tomographic concepts are able to perfectly identify no fractured areas and the main characteristics. The experiment should be repeated in the other direction, since the part of the cross-sectional area where the receiver points are located provides better results to further improve our inversion algorithm. This uneven sensitivity is a known attribute of tomographic methods, as it was shown for continuous reconstruction in [20].

The presented procedure for the estimation of variance includes noise variance and model errors. The estimated variance is usually higher than the variance of the noise, since the estimated variance includes not only the variance of the noise, but also the error of the conceptual model, e.g., the geometrical constraints (only two different fracture inclinations, fracture length distribution), the limited number of iterations, and errors of the forward model. Therefore, this method could also be applied to determine the errors of the conceptual model, if observation data without noise or with very low noise variances are used for the inversion.

Variance estimation becomes even more significant when considering, for instance, heat as a tracer. In-situ temperature observations are more prone to systematic errors and to being affected by diffusion than solutes. The simulation of thermal transport in both fracture and matrix introduces an extra level of complexity into the forward model, and requires the specification of additional physical parameters. Variance estimation could be employed to reveal initial improper settings if these parameters are set incorrectly.

The aim of this study is to extend the inversion procedure from synthetic DFN cases to practical applications with field data. To prepare this, we conclude from the synthetic DFN case study to the parameters and information we need to use this inversion algorithm in practice. In addition, we checked the limitations and the conceptual errors of this approach. However, we introduce a method to deal with that issue. One conclusion is that we need a 3D forward model to expand from the two-dimensional cases to three-dimensional DFN reconstruction. Moreover, we require information regarding the fracture sets, e.g., from outcrops. We also have to use an optical or acoustical televiewer to get the inclination and the position of the fractures adjacent to the boreholes to introduce them as source and receivers for the inversion and as boundary conditions in the forward model.