A Shakas

N Linde

T Le Borgne

O Bour

Probabilistic inference of fracture-scale flow paths and aperture distribution from hydrogeophysically-monitored tracer tests

Keywords: Fractured Rock, Hydrogeophysics, Fracture Aperture Distribution, Flow Channeling, Dispersion

Fracture-scale heterogeneity plays an important role in driving dispersion, mixing and heat transfer in fractured rocks. Current approaches to characterize fracture scale flow and transport processes largely rely on indirect information based on the interpretation of tracer tests. Geophysical techniques used in parallel with tracer tests can offer time-lapse images indicative of the migration of electrically-conductive tracers away from the injection location. In this study, we present a methodology to invert time-lapse ground penetrating radar reflection monitoring data acquired during a push-pull tracer test to infer fracture-scale transport patterns and aperture distribution. We do this by using a probabilistic inversion based on a Markov chain Monte Carlo algorithm. After demonstration on a synthetic dataset, we apply the new inversion method to field data. Our main findings are that the marginal distribution of local fracture apertures is well resolved and that the field site is characterized by strong flow channeling, which is consistent with interpretations of heat tracer tests in the same injection fracture.

Introduction

Flow and solute transport in fractured rock plays a key role in a range of applications, such as groundwater remediation (NRC, 1996;[START_REF] Berkowitz | Characterizing flow and transport in fractured geological media: A review[END_REF][START_REF] Neuman | Trends, prospects and challenges in quantifying flow and transport through fractured rocks[END_REF][START_REF] Singhal | Applied Hydrogeology of Fractured Rocks[END_REF], the design of storage facilities for spent nuclear fuel [START_REF] Tsang | Hydrologic issues associated with nuclear waste repositories[END_REF] and geothermal energy production [START_REF] Bödvarsson | Injection and thermal breakthrough in fractured geothermal reservoirs[END_REF][START_REF] Geiger | Non-Fourier thermal transport in fractured geological media[END_REF]. Fracture properties such as the distribution of fracture length [START_REF] Bonnet | Scaling of fracture systems in geological media[END_REF], surface roughness (e.g., [START_REF] Poon | Surface measurement and fractal characterization of naturally fractured rocks[END_REF][START_REF] Schmittbuhl | Scaling invariance of crack surfaces[END_REF][START_REF] Neuville | Hydrothermal coupling in a self-affine rough fracture[END_REF] and fracture aperture [START_REF] Belfield | Multifractal characteristics of natural fracture apertures[END_REF][START_REF] Lanaro | A random field model for surface roughness and aperture of rock fractures[END_REF] vary over multiple spatial scales. In geological settings where the rock matrix is largely impermeable [START_REF] Brace | Permeability of crystalline rocks: New in situ measurements[END_REF], fluid flow occurs almost exclusively within fractures. The resulting flow follows preferential flow paths and it is often characterized as channelized, leading to early tracer breakthrough and heavy tailing (e.g., [START_REF] Becker | Tracer transport in fractured crystalline rock: Evidence of nondiffusive breakthrough tailing[END_REF][START_REF] Kang | Impact of velocity correlation and distribution on transport in fractured media: Field evidence and theoretical model[END_REF]. Heterogeneous flow velocities at the fracture scale, mainly arising from local aperture disparities, represent an important driver of solute dispersion [START_REF] De Dreuzy | Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (DFN)[END_REF][START_REF] Wang | Non-fickian transport through two-dimensional rough fractures: Assessment and prediction[END_REF][START_REF] Fiori | Power law breakthrough curve tailing in a fracture: The role of advection[END_REF][START_REF] Kang | Emergence of anomalous transport in stressed rough fractures[END_REF][START_REF] Tang | Contaminant transport in fractured porous media: Analytical solution for a single fracture[END_REF][START_REF] Novakowski | A field example of measuring hydrodynamic dispersion in a single fracture[END_REF][START_REF] Moreno | Some anomalous features of flow and solute transport arising from fracture aperture variability[END_REF][START_REF] Nordqvist | A variable aperture fracture network model for flow and transport in fractured rocks[END_REF][START_REF] Bodin | Solute transport in a single fracture with negligible matrix permeability: 1. fundamental mechanisms[END_REF]. Along with the presence of aperture variations, the frequency and location of contact regions (i.e., regions of zero aperture due to fracture closure or clogging due to mineralization) have a critical impact on channeling (e.g., [START_REF] Tsang | Flow channeling in heterogeneous fractured rocks[END_REF] and fluid percolation [START_REF] Broadbent | Percolation processes: I. crystals and mazes[END_REF]). Yet, the spatial distribution of aperture within fracture planes [START_REF] Oron | Flow in rock fractures: The local cubic law assumption reexamined[END_REF] is inherently difficult to measure in situ, especially far away from boreholes.

Geophysics offers techniques for non-invasive detection of fracture geometry and fracture aperture characterization. A widely used method in this regard is the ground penetrating radar (GPR) method, in which a high-frequency radar antenna emits an electromagnetic wave that propagates through the rock matrix and is scattered at fracture locations. The recorded scattered field carries information about fracture properties, namely the aperture and orientation of the fracture as well as the electrical properties of the fracture filling (e.g., [START_REF] Bradford | Ground-penetrating radar theory and application of thin-bed offset-dependent reflectivity[END_REF][START_REF] Tsoflias | Ground-penetrating-radar response to fracture-fluid salinity: Why lower frequencies are favorable for resolving salinity changes[END_REF][START_REF] Shakas | Effective modeling of ground penetrating radar in fractured media using analytic solutions for propagation, thin-bed interaction and dipolar scattering[END_REF]. In low-loss (i.e., electrically-resistive) media, the GPR propagation distances can span several tens of meters. If the contrast in electrical properties between the fracture filling and the surrounding rock matrix is strong, clear reflections can result from fractures with apertures that are several orders of magnitude smaller than the dominant GPR wavelength (e.g., Dorn et al., 2012). The GPR method has been used for various fracture-related purposes such as evaluating methodologies for storing nuclear fuel [START_REF] Olsson | Borehole radar applied to the characterization of hydraulically conductive fracture zones in crystalline rock[END_REF], rock fall hazard [START_REF] Jeannin | Multiconfiguration GPR measurements for geometric fracture characterization in limestone cliffs (Alps)[END_REF], building safety [START_REF] Leucci | Detection of fractures from GPR data: the case history of the Cathedral of Otranto[END_REF][START_REF] Sambuelli | Estimation of thin fracture aperture in a marble block by GPR sounding[END_REF] and aquifer characterization (Dorn et al., 2012). In hydrogeophysical applications, GPR complements hydrological tests and helps to better inform the interpretation of classical hydrological data (e.g., Rubin and Hubbard, 2005, Ch. 7). During saline tracer tests, for example, GPR can provide time-lapse images indicative of tracer migration [START_REF] Ramirez | Investigation of fracture flow paths using alterant geophysical tomography[END_REF][START_REF] Day-Lewis | Time-lapse imaging of salinetracer transport in fractured rock using difference-attenuation radar tomography[END_REF]Tal-ley et al., 2005;[START_REF] Tsoflias | Ground-penetrating-radar response to fracture-fluid salinity: Why lower frequencies are favorable for resolving salinity changes[END_REF][START_REF] Becker | Comparing flux-averaged and resident concentration in a fractured bedrock using ground penetrating radar[END_REF][START_REF] Dorn | Single-hole GPR reflection imaging of solute transport in a granitic aquifer[END_REF][START_REF] Shakas | Hydrogeophysical characterization of transport processes in fractured rock by combining push-pull and single-hole ground penetrating radar experiments[END_REF]Shakas et al., , 2017)). These time-lapse images provide rich information about the dynamics of flow-and-transport during the experiment, that in turn can provide constraints on subsurface properties of interest such as aperture heterogeneity. While previous studies have focused on fracture geometry and mean aperture, they have not investigated so far the possibility to use GPR imaging to infer the aperture distribution within a fracture.

Modeling GPR propagation and scattering from heterogeneous fractures is very challenging when using the Finite-Difference Time-Domain (FDTD) formulation (e.g., [START_REF] Warren | An advanced GPR modelling framework: The next generation of gprMax[END_REF]. With FDTD it is computationally extremely demanding to model realistic fractures, for example, mm or sub-mm thin fractures in a domain of several tens of meters. Instead, most fracture-related GPR studies have relied on an analytical solution, namely the thin-bed approximation, to model EM interaction with fractures [START_REF] Tsoflias | Investigating multi-polarization GPR wave transmission through thin layers: Implications for vertical fracture characterization[END_REF][START_REF] Bradford | Ground-penetrating radar theory and application of thin-bed offset-dependent reflectivity[END_REF][START_REF] Deparis | On the use of dispersive APVO GPR curves for thin-bed properties estimation: Theory and application to fracture characterization[END_REF][START_REF] Sassen | 3D polarimetric GPR coherency attributes and full-waveform inversion of transmission data for characterizing fractured rock[END_REF][START_REF] Sambuelli | Estimation of thin fracture aperture in a marble block by GPR sounding[END_REF][START_REF] Babcock | Reflection waveform inversion of ground-penetrating radar data for characterizing thin and ultrathin layers of nonaqueous phase liquid contaminants in stratified media[END_REF]. [START_REF] Shakas | Effective modeling of ground penetrating radar in fractured media using analytic solutions for propagation, thin-bed interaction and dipolar scattering[END_REF] introduced a new methodology (the so-called effective-dipole approach) to model GPR scattering from heterogeneous fractures that is inspired by a microscopic treatment of Maxwell's equations. It is capable of modeling reflections arising from fractures with spatially-varying electrical properties and aperture along the fracture plane. In a later study, Shakas and Linde (2017) demonstrated that the thin-bed approximation is only valid for a very limited family of fracture models, namely those with smooth (large Hurst exponents) and large-scale (large correlation length) aperture variations. They proposed that explicit consideration of aperture heterogeneity is necessary for GPR-based aperture estimation even when targeting a mean aperture.

In this study, we develop a methodology to invert GPR reflection monitoring data from push-pull tracer tests to infer fracture-scale aperture variations and flow paths. For this, we combine simulations of flow and transport of an electrically-conductive tracer within a single heterogeneous fracture with the associated GPR response using the effective-dipole approach by [START_REF] Shakas | Effective modeling of ground penetrating radar in fractured media using analytic solutions for propagation, thin-bed interaction and dipolar scattering[END_REF]. More specifically, we use the local cubic law to simulate fluid flow by assigning the local fracture transmissivity based on the local aperture, while the GPR-response is determined by a semi-analytical formulation that is strongly dependent on aperture, as well as electrical conductivity and permittivity. We incorporate these simulation capabilities within a probabilistic inversion framework, using a Markov chain Monte Carlo (MCMC) algorithm, to infer a fracture's geometric properties and aperture distribution. We first apply the methodology to a synthetic example, and proceed to inversion of single-hole GPR monitoring data from a push-pull tracer test experiment conducted in Ploemeur, France (Shakas et al., 2017). For both the synthetic and real data examples, we rely on the electrically-conductive and neutrally-buoyant tracer introduced by Shakas et al. (2017). By inferring a heterogenous aperture field, we investigate if aperture fields exist for which our flow and transport modeling approach (determining where the tracer is found at a given time) and the corresponding simulated GPR-response (the GPR response from a tracer 2 Methodology

Local aperture variations

Aperture variations are generally spatially-correlated along the fracture plane [START_REF] Brown | Correlation between the surfaces of natural rock joints[END_REF] and similar aperture patterns occurring at different scales are often described using the concept of self-affinity [START_REF] Mandelbrot | Self-affine fractals and fractal dimension[END_REF][START_REF] Mandelbrot | Multifractal measures, especially for the geophysicist[END_REF]. Mathematically, self-affinity ensures that aperture fields retain statistical similarity by an appropriate scaling factor denoted by the Hurst exponent H [START_REF] Voss | Fractals in nature: From characterization to simulation[END_REF]. The exponent H can vary between 0 to 1, with H = 0 implying self-similarity (aperture patterns are statistically similar at all spatial scales) and H = 1 implying smoothly varying (multivariate normal) fields. Furthermore, H is related to the fractal dimension (D), which for 2-dimensional fields is given by D = 3 -H.

Herein, we use the following correlation function [START_REF] Adler | Fractured Porous Media[END_REF] to impose spatial correlation in fracture-wall asperity:

C h (n, m) = σ 2 h exp - n l n 2 + m l m 2 H , (1)
where the characteristic length of the spatial correlation between points separated by distance {n, m} in two perpendicular directions is given by l n and l m , and the variance of the covariance function is σ 2 h .

To create an aperture field realization along the fracture plane, we proceed by assigning the variables of the correlation function and, using the algorithm described by [START_REF] Laloy | Probabilistic inference of multi-Gaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction[END_REF], we create two independent realizations of fracture-wall asperity. We then overlay these two surfaces, which have the same statistical properties, and displace them by a mean fracture aperture. Any negative values that are created in this process are assigned an aperture of zero; this allows us to represent closed fracture sections in which no flow is allowed. This process does not mimic the mechanics of fracture generation but allows us to create a fracture aperture field with the desired statistical properties and with fracture closure (zero aperture).

In all the results that follow we discretize a single fracture of 16 m × 16 m into cells of 20 cm × 20 cm (in accordance with modeling recommendations for a 100 MHz antenna, see [START_REF] Shakas | Effective modeling of ground penetrating radar in fractured media using analytic solutions for propagation, thin-bed interaction and dipolar scattering[END_REF]) leading to aperture fields consisting of 6400 cells. We use a dimensionality-reduction algorithm [START_REF] Laloy | Probabilistic inference of multi-Gaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction[END_REF][START_REF] Hunziker | Inference of multi-Gaussian relative permittivity fields by probabilistic inversion of crosshole ground-penetrating radar data[END_REF] that allows us to represent each aperture realization with 100 dimension-reduction (DR) variables and 5 global geostatistical variables. An illustration is given in Figure 1.

Flow and transport modeling

To solve flow in a heterogeneous fracture, we use the lubrication approximation (e.g., Zimmerman and Yeo, 2013;[START_REF] Brush | Fluid flow in synthetic rough-walled fractures: Navier-Stokes, Stokes, and local cubic law simulations[END_REF]. Hence, we discretize a heterogeneous fracture into a collection of cells with each cell being treated as an equivalent parallel plate. Even if aperture can vary significantly below the cell-scale, we assume that the average aperture of each parallel plate provides an adequate description of Darcy flow. Extensions of the parallel plate model that account for more complicated geometries have been proposed [START_REF] Oron | Flow in rock fractures: The local cubic law assumption reexamined[END_REF][START_REF] Lee | Assessment of the validity of Stokes and Reynolds equations for fluid flow through a rough-walled fracture with flow imaging[END_REF][START_REF] Wang | Modification of the local cubic law of fracture flow for weak inertia, tortuosity, and roughness[END_REF]. We further assume an incompressible single-phase fluid exhibiting laminar flow.

Flow in a fracture

For laminar flow, the Navier-Stokes equations simplify to the Stokes equation. Its exact solution for a pair of parallel plates, assuming no-flow and no-slip boundary conditions, is described by the fluid velocity u (m • s -1 ):

u(z) = |∇P | 2µ z 2 - a pp 2 2 . ( 2
)
In Eq. 2, the plates are oriented along the x -ŷ plane, the aperture a pp varies along the ẑ-direction, and the pressure gradient ∇P (Pa • m -1 ) is constant along the x -ŷ direction and zero along the ẑ direction.

Integrating Eq. 2 along the ẑ-direction and multiplying by the width L pp (m) of the plate gives the volumetric flux through the parallel plate,

Q = -|∇P |L pp a 3 pp 12µ ≡ -kA|∇P | µ . (3)
The right side of Eq. 3 is equivalent to the volumetric flux derived from Darcy's Law. In 1D, Darcy's Law states that the flow of a liquid with dynamic viscosity µ (Pa • s) through a porous medium with permeability k (m 2 ) and cross sectional area A (m 2 ) is directly proportional to the pressure gradient |∇P |. The applicability of Darcy's law to the parallel plate model leads to the following analogy between permeability and aperture:

k = a 2 pp 12 . (4)
The product between permeability and area gives the "cubic law" relation, kA = Lppa 3 pp 12 . We apply the cubic law locally along a fracture's plane, an approach that is commonly termed the "local cubic law" (LCL) (e.g., [START_REF] Oron | Flow in rock fractures: The local cubic law assumption reexamined[END_REF][START_REF] Brush | Fluid flow in synthetic rough-walled fractures: Navier-Stokes, Stokes, and local cubic law simulations[END_REF][START_REF] Konzuk | Evaluation of cubic law based models describing single-phase flow through a rough-walled fracture[END_REF][START_REF] Qian | Experimental study of the effect of roughness and Reynolds number on fluid flow in rough-walled single fractures: A check of local cubic law[END_REF][START_REF] Lee | Assessment of the validity of Stokes and Reynolds equations for fluid flow through a rough-walled fracture with flow imaging[END_REF][START_REF] Wang | Modification of the local cubic law of fracture flow for weak inertia, tortuosity, and roughness[END_REF]. The validity of the LCL depends on the following conditions [START_REF] Oron | Flow in rock fractures: The local cubic law assumption reexamined[END_REF]:

σ {u,d} pp a pp << 1, a pp 2L pp 2 << 1, and Re • max σ pp a pp , a pp 2L pp << 1, (5)
with σ pp being the standard deviation of the aperture along the extent of a parallel plate. The third condition in Eq. 5 includes the Reynolds number, defined as Re ∼ = ρU app 2µ

with U (m • s -1 ) being the characteristic velocity of the fluid flowing through the parallel plate. In our simulations we define L pp = 20 cm, but we do not directly assign values for σ pp and a pp ; instead, we create the aperture field directly, and hence all local apertures (a pp ) along the fracture plane, with the algorithm described in sec. 2.1.

This procedure results in the aperture distribution of a pp seen in Fig. 6(c). The validity of the second inequality can be verified by using the maximum value of the posterior distribution for a pp , which is 8 mm (introduced later in Fig. 6). This results in 0.0004<<1. To check the third inequality we need a characteristic flow velocity in order to compute the Reynolds number. For this, we consider the 3 mm•s -1 average flow velocity estimated by [START_REF] Shakas | Hydrogeophysical characterization of transport processes in fractured rock by combining push-pull and single-hole ground penetrating radar experiments[END_REF] for an almost identical experiment. This results in Re = 4.29 and the limiting factor is, thus σ pp /a pp << 0.23. We see that the LCL is only valid when fracture aperture heterogeneity is negligible below L pp .

Solute transport

Once the fluid velocity u has been determined at each location along the fracture plane, it can be used to compute the transport of a solute by the advection-diffusion equation:

∂ ∂t (c) + ∇ • [cu -D m ∇c] = 0. ( 6
)
The concentration of the solute in Eq. 6 is denoted by c (kg • m -3 ) and D m is the molecular diffusion coefficient for which we use a fixed value of D m = 2 × 10 -9 m 2 • s -1 [START_REF] Holz | Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate 1H NMR PFG measurements[END_REF]. We ignore Taylor dispersion in the following and assume that the impact of numerical dispersion is negligible. Since our GPR response is mainly sensitive to regions of high salinity, we expect that dispersion has only a minor effect in the present study (see Fig. 4).

Flow and transport in a rough fracture

For modeling flow and transport of a tracer in a single heterogeneous fracture, we modified the MaFloT-2D code that was developed for porous media by [START_REF] Künze | MaFloT-matlab flow and transport[END_REF]. MaFloT-2D is a finite volume algorithm that simulates density-driven flow, hydrodynamic dispersion, diffusion and advection in a 2D cross section of a porous medium by assigning the porosity and permeability of each cell in the 2D section.

It solves the mass balance equation for a single incompressible fluid completely saturating a porous matrix to obtain the pressure distribution. Using the pressure solution, it computes the Darcy velocity of the fluid within the plane and use it to compute the concentration at each time step by discretizing the advection-diffusion equation in both time and space.

Our modifications to MaFloT-2D are the following:

1. Assign the permeability of each cell using the permeability-aperture relation, k = a 2 pp 12 .

2. Replace the porosity of each cell by the local aperture a pp .

3. Use the minimum value as the interface area between cells, A 1,2 min = L pp × min{a 1 pp , a 2 pp } where the superscripts 1 and 2 denote two adjacent cells.

4. Introduce contact cells (zero aperture) as no-flow zones.

We assume flow and transport to be advection-dominated and assign either Dirichlet (zero-pressure) or

Neumann (no-flow) boundary conditions at the outer boundaries as shown in Figure 2. More specifically, we assign no-flow to the bottom, left and right sides of the main fracture and zero pressure on the top, thus forcing the tracer to migrate upwards. To simulate the injection and withdrawal of a tracer during a push-pull test, we assign dynamic inflow and outflow boundary conditions at the injection location, which we always place at the center-bottom of a fracture (corresponding to the borehole location). This was done to obtain a similar flow configuration as the one observed in our field-experimental results that demonstrate upward movement of tracer from the injection depth [START_REF] Shakas | Hydrogeophysical characterization of transport processes in fractured rock by combining push-pull and single-hole ground penetrating radar experiments[END_REF](Shakas et al., , 2017)).

Ground penetrating radar modeling

To model the propagation of the GPR source wavelet in the rock matrix and the subsequent interaction with the tracer-filled fracture, we use the effective-dipole methodology introduced by [START_REF] Shakas | Effective modeling of ground penetrating radar in fractured media using analytic solutions for propagation, thin-bed interaction and dipolar scattering[END_REF]. We allow for spatial variations in aperture and fracture-filling properties following Shakas and Linde (2017) and use the simulated distribution of tracer obtained from the MaFloT-2D simulation to compute time-lapse changes in the GPR signal. We do this at each measurement location by computing the difference between the later traces (when tracer is already in the fracture) with the background trace (before the tracer experiment starts); a visualization of this process is given by Fig. 4 in Dorn et al. (2012).

We vary the shape of the GPR source spectrum using the four parameters of the generalized Gamma distribution. Given that the actual GPR source energy is unknown in field situations, we normalize each forward GPR simulation such that the total reflected energy in every GPR simulation matches the total reflected energy in the observed data. We compute the total energy by summing the absolute values of the complex-valued GPR traces expressed in the frequency-domain; the latter are obtained directly from the simulations (the forward model is constructed in the frequency domain) and through a Fast Fourier

Transform (f f t) of the real data. For a single trace, the total energy

(T E) is T E = N 1=1 |f f t(t i )| = N 1=1 (x 2 i + iy 2 i ) 0.5
, where t i is the time-varying (real-valued) amplitude of the electric field for a trace of N samples. The total reflected energy is simply the sum of the total energy over all the traces considered in a given simulation.

Tracer and rock properties

Throughout this work we assume the tracer to be wethanalt, a neutrally-buoyant and electricallyconductive tracer introduced by Shakas et al. (2017). Wethanalt is a mixture of saline tracer and ethanol; the latter allows to adjust the buoyancy of the saline tracer to match the ambient formation water, while the former ensures a strong electrical conductivity contrast with the ambient water. Further information about how to prepare wethanalt and experimental findings from its use are presented in Shakas et al. (2017). We define the following fluid properties for each parallel plate:

• Electrical conductivity σ (S • m -1
), which we link to salt concentration c (g • L -1 ) using Eq. 9 of [START_REF] Sen | Influence of temperature on electrical conductivity on shaly sands[END_REF].

• Dynamic viscosity µ (mPa • s) which varies linearly between the value for ambient water (µ = 1) to the value for wethanalt (for our experiments the wethanalt solution has µ = 2.8). This approximate linear relation between water and ethanol over this range is supported by controlled lab experiments presented by [START_REF] Hammond | CRC Handbook of Chemistry and Physics[END_REF].

• Relative electrical permittivity ε r , vary according to the water/ethanol ratio that is computed from the tracer concentration in each plate. The values range linearly from wethanalt (ε r = 53) to ambient water (ε r = 79) following [START_REF] Sato | Dielectric relaxation processes in ethanol/water mixtures[END_REF].

Note that all the parameters introduced above are temperature dependent and are here assumed at their value at 20 o C (the ambient water temperature at the field site is 16 o C). Furthermore, we do not consider density effects since the density of wethanalt is by design the same as the ambient water. In addition to the electrical properties of the fracture filling, also the electrical properties of the rock matrix need to be defined. In both the synthetic and real data inversions, we set the same realistic prior ranges for these parameters that were obtained from existing studies of the rock properties at the Ploemeur field site [START_REF] Belghoul | Caractérisation pétrophysique et hydrodynamique du socle cristallin[END_REF] where we conducted the experiments [START_REF] Shakas | Hydrogeophysical characterization of transport processes in fractured rock by combining push-pull and single-hole ground penetrating radar experiments[END_REF](Shakas et al., , 2017)).

The prior parameter ranges for the parameters of interest are presented in Table 1. All the listed parameters have a bounded uniform prior distribution, except for the mean aperture α and standard deviation σ h for which we used a bounded log-uniform prior to account for the fact that these parameters may vary over several orders of magnitude (sec. 2.1). The prior on each DR variable is given by an uncorrelated standard normal distribution.

Model geometry

We consider the same model geometry for both the real and synthetic data inversions (Figure 2). Our previous experience with the considered field data has shown that the tracer moves upwards after injection [START_REF] Shakas | Hydrogeophysical characterization of transport processes in fractured rock by combining push-pull and single-hole ground penetrating radar experiments[END_REF](Shakas et al., , 2017)). Nevertheless, televiewer data indicates that the injection fracture is oriented sub-horizontally with a deviation of 15 o from the horizontal plane (fracture B1-4 in Le Borgne et al.

(2007)). These observations suggest that the sub-horizontal fracture, which is not seen in the GPR reflection data, is connected to a neighboring sub-vertical fracture. To account for this, we define an initial sub-horizontal fracture (dip = 15 o ) in which we inject the tracer. This fracture folds into a subvertical fracture, whose dip and azimuth we infer during inversion, together with the folding distance (F ) from the injection location that can range from 20 cm to 4 m (in steps of 20 cm given the considered cell size). The fixed fracture side length and width (both 16 m) allow for enough fracture area to reproduce the upward movement (approximately 10 m) of the tracer observed in the experiments. We do not compute the GPR reflections arising from the sub-horizontal fracture as they do not present a significant response.

During inversion, we only accept models with an aperture distribution that is hydraulically connected from the injection location to the top of the fracture. We assign fixed borehole trajectories, for both the injection and monitoring boreholes, by utilizing previous borehole deviation logging and cross-hole GPR tests (Dorn et al., 2012). This allows us to define the position and orientation of the GPR transmitter and receiver at each simulation step.

Probabilistic inference

Our primary focus is to infer the geometric properties of the sub-vertically-oriented fracture, namely its orientation and aperture distribution, as well as the folding distance F . To do this, we rely on the DREAM (ZS) algorithm [START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing[END_REF]; a Markov chain Monte Carlo (MCMC) algorithm that uses differential evolution updating of the model parameters to efficiently sample their posterior probability densities. We use the DREAM (ZS) algorithm in parallel over 24 processors and assign one chain for each processor.

As data, we use the absolute value of the complex-valued coefficients of each frequency component of the Fourier transformed GPR differences (see section 2.3). The use of absolute values results in smooth data that are easier to fit, compared with highly oscillating waveforms, but without information about the phase of the incoming electric field. However, the absolute values are directly affected by the spatial distribution of tracer in the fracture and are informative of the aperture distribution along the fracture plane. We assume that each real and imaginary part has normally-distributed noise with zero mean and standard deviation σ n . The resulting likelihood function (L) characterizing the absolute values is a product of Gamma distributions (see derivation in Appendix A). The maximum-likelihood fit when using the gamma distribution results in a model that has an average misfit σ n between the observed and modeled data. In contrast, when using the normal distribution the maximum-likelihood model favors exact overlap between observed and modeled data (zero-misfit). The likelihood function is:

L = M m=1 γ(∆d m , σ n ) = M m=1 ∆d m σ 2 n exp - ∆d m σ n 2 . ( 7
)
where ∆d m is the difference in absolute values between the mth observed (d m o ) and simulated (d m s ) data (in vector notation:

d o = {d 1 o , d 2 o , .., d M o }).
log L = M m=1 log γ(∆d m , σ n ) = M m=1 -2log(σ n ) + log(∆d m ) - ∆d m σ n 2 . ( 8
)

Parameters of interest

To infer the fracture aperture distribution, we invert for the mean aperture α, standard deviation σ h , two correlation lengths l n and l m and Hurst exponent H, as well as 100 DR variables (sec. 2.1). For the fracture geometry, we invert for the dip θ and azimuth ϕ as well as the folding distance F (sec. 2.5). For the electric properties of the rock, we invert for both electrical permittivity ε r and conductivity σ (sec.

2.4) and for the four source parameters [START_REF] Shakas | Effective modeling of ground penetrating radar in fractured media using analytic solutions for propagation, thin-bed interaction and dipolar scattering[END_REF]. Additionally, we apply a Hierarchical Bayes scheme and consider the noise level (σ n in Eq. 8) as a hyper-parameter that varies in the range

σ min ≤ σ n ≤ 10 • σ min .
Here, σ min is the minimum level of noise expected in the data; this corresponds to the known added noise for the synthetic case and it is computed for the real data by examining regions in the GPR difference sections where no changes are expected to occur (i.e., where there are no reflections).

We assign a low log-likelihood (by replacing ∆d m with ∆d m = 10 • d m o in Eq. 8) to all proposed aperture models for which the injection occurs within a zero aperture region, or if there is no percolation from the injection cell to the fracture's open boundary. In the following section we present the inversion results obtained for the synthetic and field-based examples.

Results

We begin this section by presenting inversion results from the synthetic test case before proceeding to those of the field experiment. In the synthetic example, we chose the experimental parameters, both for the hydrologic test and the geophysical monitoring, to match those used in the actual field experiment.

For the latter we consider the push-pull tracer test (a) from Shakas et al. (2017); further information about combined push-pull and GPR experiments at this site can also be found in [START_REF] Shakas | Hydrogeophysical characterization of transport processes in fractured rock by combining push-pull and single-hole ground penetrating radar experiments[END_REF].

We assign a constant pumping rate of 2.7 L • min -1 during the injection, chasing and withdrawal phases and model the push-pull test by injecting wethanalt with electrical conductivity of 35 mS • cm -1 , corresponding to 44 g • kg -1 of salt, into a fracture filled with fresh water of 0.7 mS • cm -1 for 33 min, followed by a chasing with formation water for another 33 min before reversing the flow.

We use GPR data from 43 locations with mid-points that are 50 cm apart along the borehole; the actual field acquisition was made at 5 cm intervals but we consider only every 10th trace to gain computational time. We choose 6 difference sections from the field dataset that showed strong changes in the GPR reflections. The times at which these six sections were measured are {t 1 : 23 to 26, t 2 : 29 to 32, t 3 : 37 to 40, t 4 : 44 to 47 t 5 : 51 to 54 and t 6 : 59 to 62} minutes, measured from the start of the tracer injection.

Here we present tracer simulation snapshots and thus refer to three time instants as t 2 = 30, t 4 = 46 and t 6 = 61, that we consider representative of sections t 2 , t 4 and t 6 respectively. In practice, during the simulations we compute the GPR traces using the tracer distribution snapshot that corresponds to the instant at which the antenna system is present at a given location along the borehole. Nevertheless, these tracer distribution snapshots do not vary greatly during the acquisition of a complete GPR section.

Additionally, from each FFT-transformed trace, we select samples at 80 linearly spaced frequencies in the range 0 ≤ f ≤ 190 MHz. The final dataset for both field and synthetic data consists of 20640 (43 traces, 80 samples, 6 sections) datapoints. For all the simulation results presented herein, the origin of the coordinate system (x = 0, y = 0, z = 0) is the tracer injection point. In the field experiment this corresponds to a depth of 77.8 m.

Inversion results for synthetic data

For the synthetic test case we create a fracture realization which we refer to as the true model. The geometric properties of the true model are presented in Table 1 and its geometry is shown in Fig. In Fig. 6(b) we plot the tracer probability map after inversion. This map is computed as follows: (1) the tracer simulation is run for each of the last models of the 24 MCMC chains;

(2) we assign a value of 1 to each cell of the fracture for which the tracer concentration is above 4.4 g • L -1 (10 % of the injected salt

Inversion results for field data

We now present the inversion results obtained for the field data. The prior parameter ranges and pseudoposterior statistics are presented in Table 1; as mentioned above, all settings are the same as for the synthetic example. Inversion of the field data leads to three model categories that are represented in Fig. 8. A planar view of these models along with the corresponding tracer snapshots at three time instances t 2 , t 4 and t 6 are shown in Fig. 9. The inferred models favor upwards flow of the tracer with some latency:

the tracer tends to stay between 8 and 12 m at t 2 and t 4 (second and third row, Fig. 9) and some mass remains at the end of the chasing period. Additionally, the first two types of inferred models show very similar structure, but they are mirrored along the azimuth; this results from the inability of classical omni-directional borehole GPR antennas to delineate the azimuthal direction of the reflections [START_REF] Olsson | Borehole radar applied to the characterization of hydraulically conductive fracture zones in crystalline rock[END_REF]. The third model differs from the first two in that it has a smaller correlation length along the fracture width (l n ), namely 2.7 m instead of roughly 6 m for the other two models, but it still leads to a channelized flow path.

In Figs. 10(a,d,g) we plot the observed GPR difference sections corresponding to time instances t 2 , t 4 and t 6 . The corresponding GPR simulations for the ML model are shown in Figs. 10(b, e, h), and the absolute differences between the field data and the simulated ML data are shown in Figs. 10(c,f,i).

The ML model reproduces the main apects of the data (vertical position and spectral content), while the details are not that well reproduced as for the synthetic example. In Fig. 10(j) we plot the evolution of the likelihoods for the 24 chains during the MCMC inversion, as well as the maximum and minimum error lines. Compared to the synthetic inversion, the ML model barely reaches the maximum error bound.

Also, the log-likelihood is still slowly increasing at the end of the MCMC chains and the final 24 models can consequently only be considered as an approximate pseudo-posterior.

In Fig. 11(a) we plot the tracer probability map computed using the last 24 MCMC model (see explanation in sec. 3.1). Alongside in Fig. 11(b) we plot the prior and pseudo-posteriors pdf 's of the local aperture with the pseudo-posterior being focused towards smaller local aperture values. Figure 12 shows the prior and pseudo-posterior histograms of the parameters tabulated in Table 1. As in the synthetic inversion, the pseudo-posteriors show more peaked distributions than the priors, suggesting that the data are informative in constraining these parameters.

Discussion

The principal aim of this study is to infer the aperture distribution and assessing its impact on tracer transport. For instance, does the GPR data allow us to distinguish between an aperture distribution for which the tracer spreads radially and remains close to the injection depth and a fracture exhibiting channeling in which the tracer quickly migrates upwards? The prior distribution does not allow us to distinguish between these two cases.

In this section we compare the results obtained by the synthetic and field-data inversions. In doing so, our aim is to assess the field-data inversion, for which significant modeling and experimental uncertainties exist, in light of the synthetic inversion, for which the same forward model is used to both create the data and infer the fracture properties, and for which there is no experimental uncertainty except for Gaussian noise added on the GPR data. By modeling uncertainties we refer to inadequacies of the forward model to capture the real phenomena under study (e.g., simplified fracture geometry, local cubic law, ignoring dispersion and matrix diffusion) and by experimental uncertainties we refer to data errors, sensor position uncertainties, and in general everything related to the measurement process that can not be exactly monitored and reproduced in the simulations.

A major challenge in this work is the non-linearity of the forward problem. We use the difference between two covariance functions to create a fracture aperture distribution. Flow and transport are then solved numerically to simulate the tracer distribution in the resulting fracture. Using the tracer distribution and the fracture geometrical properties, we then compute GPR reflection data and use these data in a MCMC inversion framework to infer the fracture properties. It is evident that this sequence of non-linear modeling steps leads to a highly non-linear inference problem with multiple interacting parameters.

Synthetic inversion

For the synthetic example, all inferred fracture models result in tracer that moves upwards in the center of the fracture (last row in Fig. 4) which is also the case for the true model. Also, in all models the tracer forms a semi-circular "ring" around the injection location, with parts of the tracer remaining at the same height as the injection location even at the end of the chasing period. This form of tracer flow is described as radial, that is, the tracer spreads in a radius around the injection location. Radial flow suggests that aperture pathways are smoothly distributed around the injection location, and no preferential direction exists.

The good agreement between the true and inferred tracer distributions is also reflected in the GPR difference sections (Fig. 5): both the true and ML models show an initial strong reflectivity change ) that later moves upwards (Figs. 5(d-e)) and strongly diminishes at the last time step ). The final GPR images for both the synthetic data and the ML model show reflected energy that stays at low altitudes and frequencies (around 7 m and 70 MHz), originating from the tracer that remains at similar depths as the injection location. Overall, the radial flow pattern of the tracer is seen in the GPR reflection data as a smooth amplitude ellipse that spreads in height and decreases in amplitude during the chasing period.

Field data inversion

For the field-data inversion, we observe that the final models (Fig. 8) favor an aperture distribution that causes the tracer to be channeled away from the injection point. In Fig. 9 we note that the tracer initially (second row) either occupies the lower left or the lower right part of the fracture plane, but is subsequently pushed (third row) such that, at the last time instance shown (fourth row), there is not much tracer remaining close to the injection location. This is also seen in the GPR reflection images presented in Fig. 10, in which the field-data (first column) shows reflections at the end of the tracer injection that span a large depth range (from 0 m to roughly 20 m in height) that are later pushed during chasing to the upper region (Fig. 10(d)) and finally end up only occupying the upper region (Fig. 10(g)).

This pattern is also reproduced by the ML model (second column in Fig. 10).

In order to obtain a model that reproduces the patterns observed in the field-based GPR reflection data (Fig. 10), the inversion favors a fracture aperture distribution that allows the complete upwards migration of the tracer during the pushing phase of the experiment. This is not the case for the synthetic inversion, for which the tracer spreads in a more radial manner (Fig. 4). In the field-data, the GPR reflections The ability of the field-data inversion to infer a strong degree of flow channeling in the fracture plane is an important finding owing to the role of flow channeling in controlling dispersion and mixing processes in fractures. This is consistent with indirect evidence obtained from previous results of heat tracer tests at the site, which suggested that heat recovery was characteristic of a highly channeled flow [START_REF] Klepikova | Passive temperature tomography experiments to characterize transmissivity and connectivity of preferential flow paths in fractured media[END_REF][START_REF] Klepikova | Heat as a tracer for understanding transport processes in fractured media: Theory and field assessment from multiscale thermal pushpull tracer tests[END_REF]. Flow channeling not only affects the diffusion of heat from the fracture to the matrix, but it also has important consequence for solute transport by inducing both fast transfer times through preferential flows and very long residence times in the remaining slow-flow areas of the fracture.

Aperture distribution

Another important parameter that can be constrained from our inversion results is the local aperture distribution. In the synthetic inversion, we found that the local aperture distribution is well retrieved (Fig. 6(c)). Unlike the prior, the posterior distribution reproduces the peak, variance and skewness of the true distribution. The corresponding pseudo-posterior local aperture distribution obtained from the fielddata inversion (Fig. 11(b)) favors smaller apertures and has a shape that is strikingly different than the prior distribution. Specifically, the pseudo-posterior suggests that the local apertures are distributed over a much narrower range, and that local apertures above 5 mm are highly unlikely. The tracer probability map shown in Fig. 11(a) suggests that the tracer is equally likely to follow a left or right path (but not an direct upwards movement from the injection location); this reflects the fact that single-hole GPR reflection data is insensitive to the azimuth of the reflection (see [START_REF] Olsson | Borehole radar applied to the characterization of hydraulically conductive fracture zones in crystalline rock[END_REF] for more information).

In addition to the two modes of the posterior distribution indicating an initial sideways motion, we also find that the tracer has to move through a bottleneck at roughly 5.5 m above the injection location.

Reliability of the pseudo-posterior probability distributions

The mean and standard deviation of the pseudo-posterior probability distributions are presented in Table 1. For the synthetic data inversion (see Fig. 7), we note that the pseudo-posteriors overlap with the true parameter valuse. Furthermore, the pseudo-posteriors are more peaked, suggesting an information gain compared to the priors. For the field-data inversion (Fig. 12), we note that all pseudo-posteriors are considerably more narrow than their priors. The pseudo-posteriors suggest that the main fracture has a positive dip (Fig. 12(c)), leaning towards the GPR monitoring borehole. This fracture is likely to be connected to the sub-horizontal fracture at approximately 3 m away from the injection location. Also, the two integral scales (Fig. 12 (f) and (g)) are of similar magnitude. The Hurst exponent (Fig. subfigure 12(e)) is likely to be larger than ∼ 0.8, suggesting aperture patterns that are more similar to a multi-Gaussian distribution. For fracture surfaces, Hurst exponents of this magnitude have been supported by experimental evidence [START_REF] Brown | Correlation between the surfaces of natural rock joints[END_REF]. Unfortunately, computational resources did not allow for a formal convergence of the posterior distributions. The best approach to achieve convergence in future studies is probably offered by parallel tempering (see [START_REF] Laloy | Merging parallel tempering with sequential geostatistical resampling for improved posterior exploration of high-dimensional subsurface categorical fields[END_REF] for a demonstration of the significant improvements offered by parallel tempering in data rich environments). Basically, when using large data sets (20640 datapoints in this case) with high signal-to-noise-ratios, it becomes highly unlikely that a traditional MCMC algorithm accept uphill steps (models with lower likelihood), which implies that there is a large risk to remain stuck in local minima. Parallel tempering [START_REF] Earl | Parallel tempering: Theory, applications, and new perspectives[END_REF] circumvents this problem by allowing for information exchange with other MCMC chains that move more freely through the model space.

These posterior estimates are also affected by strong modeling assumptions; the main ones are summarized below. In terms of GPR modeling, the main assumption affecting our results is probably that we do not account for fracture topography (we account for aperture variations, but assume that the fracture is planar). This could be accounted for by using appropriate meshing techniques and including dipole coupling in the GPR modeling (Shakas and Linde, 2017). The most important assumption is probably related to the fracture geometries considered. The two intersecting fractures are likely to have a much more complex intersection than what is assumed here and additional fractures might play an important role. [START_REF] Dorn | Conditioning of stochastic 3-D fracture networks to hydrological and geophysical data[END_REF] considered the integration of GPR and hydrological data in the context of discrete fracture networks, but doing this was outside the scope of the present work. Furthermore, we ignore matrix diffusion, hydrodynamic dispersion, macro dispersion due to unaccounted heterogeneity below the discretization scale (20 cm) and we assume that the cubic law is valid locally. We expect that the impact of these assumptions are small in the present study compared with the much larger errors caused by the simplified fracture geometry model.

Conclusions

This study investigates how single-hole GPR reflection monitoring of push-pull tracer tests may provide information on fracture scale transport pathways and aperture distribution. Such properties, which can generally only be inferred indirectly (and approximately) by tracer test interpretation, play a key role in driving solute dispersion and heat transfer in fractured rocks. To target this challenge, we established a novel inverse modeling framework combining fluid flow, transport and electromagnetic wave solvers at the fracture scale with a Markov chain Monte Carlo algorithm. After demonstrating the performance of this framework for a synthetic case study, we apply it to a field experiment that used an electrically-conductive and neutrally-buoyant tracer.

From the synthetic test case we demonstrate that distinctively different aperture distributions can be found that are in strong agreement with the data. We also find that the large data volumes considered and the strong non-linearity of the involved forward solvers prohibit efficient mixing of the MCMC chains.

With the available computational budget, we were unable to sample the full posterior distribution; instead, we rather sample a pseudo-posterior that provides insights about possible fracture configurations. All the final models capture the same main attributes of the tracer migration, which for the synthetic case is characterized as radial flow. On the contrary, for the field-data inversion our final models favor a channelized tracer migration, which is in agreement with heat-tracer tests performed on the same fracture.

To the best of our knowledge, this is the first time that geophysical data has been used to infer aperture heterogeneity patterns at the fracture scale, and to provide information about their impact on flow and transport at this scale. Our main findings are that we are able to infer the marginal distribution of local apertures and distinguish between radial and channelized transport at the fracture scale, two regimes that have fundamentally different dynamics in terms of dispersion, mixing and contact area between solute and rock.

Our GPR forward model can simulate reflections arising from multiple fractures, while the flow and transport simulator used is limited to one fracture. Here, we simulate in a simplified way the interaction of two fractures by assigning two dip angles to a single fracture plane. In the future, it would be important to couple the GPR forward model with a flow and transport solver that is capable of modeling the tracer migration within a fracture network. This would allow for inference of properties characterizing fracture networks, and provide further insights to site-specific processes that are relevant, for instance, to groundwater remediation and geothermal heat extraction. For future investigations, it would be interesting to test the impact of local dispersion, which we neglected in this study. While for the considered test we do not expect local dispersion to have a significant impact on the results, for larger tracer transport distances it would certainly have an impact.

likelihood value (ML), obtained by setting the derivative with respect to w to zero, gives: ∂g(w) ∂w = 0 → ML(w) = σ n (13)

That is, the likelihood function in Eq. 12 favors an average misfit between simulated and observed data that is equal to the standard deviation of the noise in both the real and imaginary values of the transformed GPR trace.

Figure 1 :

 1 Figure 1: Fracture aperture heterogeneity is created by subtracting, (a) an independent (16 m × 16 m) spatial random field describing fracture-wall asperity, from (b) another independent realization with the same geostatistical properties. These fracture-wall asperity fields are drawn from a geostatistical model with zero mean, σ h = 0.5 mm, H = 0.8, l n = 3 m (horizontal) and l m = 5 m (vertical). (c) The resulting aperture field realization is obtained after a mean aperture of 2 mm (this example) has been added after the subtraction and any resulting intersecting regions have been set to zero aperture.

FFigure 2 :

 2 Figure 2: Schematic representation of the fracture model used in both the synthetic and real data examples considered. Injection (green disk) occurs in a sub-horizontal fracture (black) that intersects the injection borehole (light blue). The sub-horizontal fracture has variable length (F ) that we invert for. The tracer flows into a sub-vertical fracture whose dip (θ) and azimuth (ϕ) we also invert for. The fracture has no-flow boundaries (yellow) except at the top (dark blue) where zero-pressure is applied. GPR reflection data are collected in a nearby borehole (red) located ∼ 6 m from the injection borehole. Precise borehole trajectories inferred from borehole deviation logging and cross-hole tests are considered in the modeling. This simplified fracture description allows us to explain the main large-scale behavior.

  Both d m s and d m o consist of real and imaginary components, where d m = x m + iy m and |d m | = (x m ) 2 + (y m ) 2 . When using MCMC, it is numerically beneficial to evaluate model proposals by using the log-likelihood:

3

  (a), while a plane view of the fracture is shown in Fig.4(a). The parameters were chosen to be representative of values measured in the field (see http://hplus.ore.fr/en/ploemeur/publications-ploemeur for a comprehensive list of publications).The simulated tracer distributions for the three time instances t 2 , t 4 and t 6 are shown inFigs. 4(e, i, m).

Figure 3 :

 3 Figure 3: (a) The geometry of the true model and (b, c, d) the three main categories of models obtained after inversion of the noise-contaminated synthetic data, with (b) the maximum likelihood (ML) model. The aperture variations along the fracture plane are shown and the red line indicates the orientation of the GPR monitoring borehole. The injection location is defined as the origin of the coordinate system.

Figure 4 :

 4 Figure 4: Fracture aperture realizations with corresponding snapshots of tracer distributions. The first row shows the aperture distribution along the fracture plane for the (a) true model and the three categories of models (b), (c) and (d) inferred by the inversion, with (b) being the maximum likelihood model. Below each model, the tracer distribution along the fracture plane is shown at time instants t 2 , t 4 and t 6 . The injection location is indicated with a (filled green) circle and the semi-transparent area indicates the sub-horizontal part of the fracture, from which no GPR reflections are computed.

t 2 ,Figure 5 :

 25 Figure 5: Absolute amplitude of the GPR difference sections as a function of frequency and height from the injection location: (a, d, g) the true noise-contaminated model response, (b,e, h) the maximum likelihood (ML) response and (c, f, i) their differences at times t 2 , t 4 and t 6 , respectively. Each row can be compared to its representative tracer snapshot in Fig.4.(j) The evolution of the log-likelihood of the 24 MCMC chains; the dotted lines indicate the log-likelihood corresponding to the upper and lower data error bounds.

Figure 6 :

 6 Figure 6: (a) The region of the true model that is occupied by tracer at some point during the considered simulation period. (b) Corresponding probability map given as an average over the 24 final models of each MCMC chain. (c) The true, prior and pseudo-posterior distributions of the local fracture aperture (a pp ).

Figure 7 :

 7 Figure 7: Histograms of the pseudo-posteriors for the fracture parameters considered in the inversion of the synthetic data, namely, (a) mean aperture α, (b) standard deviation of aperture σ h , (c) dip θ and (d) azimuth ϕ, (e) Hurst exponent H, correlation length along the fracture (f) width l n and (g) length l m , and (h) the fold distance from injection F (see sec. 2.1 for more details). The true parameter values are shown along the horizontal axis with a (green) dot, the prior is shown with a solid (red) line and the (blue) pseudo-posterior histograms are obtained from the final 24 MCMC realizations.

Figure 8 :

 8 Figure 8: The model geometry for the three (a, b, c) main types of models obtained after MCMC inversion of the field data, with (a) being the ML model. The aperture variations along the fracture plane are shown and the red line indicates the orientation of the GPR monitoring borehole. The injection location is defined as the origin of the coordinate system.

Figure 9 :

 9 Figure 9: (a, b, c) The first row shows the aperture distribution along the fracture plane for the three types of models obtained by the MCMC inversion with (a) being the ML model. For each model, its corresponding column shows the simulated tracer distribution along the fracture plane at time instances t 2 , t 4 and t 6 . The injection location is indicated with a (filled green) circle and the semi-transparent area indicates the sub-horizontal part of the fracture, where no GPR reflections are computed.

Figure 10 :

 10 Figure 10: Absolute amplitude of the GPR difference sections as a function of frequency and height from the injection location:(a, d, g) the field-data, (b, e, h) the maximum likelihood (ML) response and (c, f, i) their differences at times t 2 , t 4 and t 6 , respectively. Each row can be compared to its representative tracer snapshot in Fig.9. (j) The evolution of the log-likelihood of the 24 MCMC chains; the dotted lines indicate the log-likelihood corresponding to the upper and lower data error bounds.

Figure 11 :

 11 Figure 11: (a) Probability that a given region in the model is occupied by the tracer at some point during the considered simulation period as given by an average over the 24 final models of each MCMC chain. (b) the prior and pseudo-posterior distributions of the local fracture aperture (a pp ).

Figure 12 :

 12 Figure 12: Pseudo-posterior histograms for the inferred fracture generation parameters used in the inversion of the field data, namely, (a) mean aperture α, (b) standard deviation of aperture σ h , (c) dip θ and (d) azimuth ϕ, (e) Hurst exponent H, correlation length along the fracture (f) width l n and (g) length l m , and (h) the fold distance from injection F (see sec. 2.1 for more details). The prior is shown with a solid (red) line and the pseudo-posterior histograms are created using the last parameter values in each MCMC chain.

(

  Fig. 10, first column) are initially weak (Fig. 10(a)), and strengthen as they move upwards from the injection location (Figs. 10(d) and (g)) as the chasing continues. This pattern is reproduced in the ML model (Fig. 10(a), second column) as well as in the other model families. In contrast to the radial tracer pattern observed for the synthetic data, the observed GPR reflection attributes can only be reproduced by a tracer that exhibits channelized upwards flow.

Table 1 :

 1 Table with results from the synthetic (subscript s) and field (subscript f ) data inversions. The prior ranges are the same for both inversions, and the pseudo-posterior mean and standard deviation are shown for the synthetic and field data. The prior type is either uniform (u) or log-uniform (l-u). For the synthetic case, the true parameters are also shown. The parameters that define the fracture geometry consist of the dip θ, azimuth ϕ, mean aperture α, standard deviation σ h , the Hurst exponent H, two characteristic correlation lengths l n and l m , as well as the folding distance F . Additionally, the electrical permittivity ε r and conductivity σ of the rock matrix are considered unknown. Results for the 100 DR variables are not shown.

	Parameter	Prior	Prior	True s	mean s	sd s	mean f	sd f
			Type
	θ ( o )	75,105	u	90	92.01	1.47	101.45	2.04
	ϕ ( o )	-22.5,22.5	u	0	3.05	8.16	1.26	18.09
	α (mm)	0.1,10	l-u	2	1.36	0.55	0.21	0.06
	σ h (mm)	0.03,3	l-u	0.15	0.16	0.03	0.14	0.07
	H (-)	0.5,1	u	0.8	0.93	0.07	0.93	0.09
	l n (m)	1,10	u	2	1.40	0.30	3.62	1.35
	l m (m)	1,10	u	6	6.30	1.95	3.66	1.09
	F (m)	0.2,4	u	2	2.33	0.32	3.04	1.03
	ε r (-)	5.3,7.4	u	7	5.77	0.23	6.68	0.51
	σ (µS • m -1 )	3,300	u	100	221.38	118.81	177.47	113.98

Acknowledgments

This research was supported by the Swiss National Science Foundation under grant 200021-146602 and by the French National Observatory H+ (hplus.ore.fr/en). The data and code required to reproduce the results are available from the first author upon request and the tracer test data are available from the H+ database.

Appendices A Likelihood function derivation

Here, we derive the likelihood function that appears in Eq. 7. To begin, we assume that the real and imaginary value pairs for each frequency component are characterized by zero-mean and normally distributed uncorrelated noise with standard deviation σ n . We denote the respective noise distributions with X = N (0, σ n ) and Y = N (0, σ n ) for the real and imaginary parts, respectively.

The absolute value of a complex number z, where z = x + i y, is |z| = x 2 + y 2 . Since we are interested in the error distribution of the absolute values, we need to propagate the errors in X and Y to the error of the absolute values Z.

The chi-squared (χ 2 ) distribution describes the sum of squares of k normally-distributed variables:

where Γ(n) ≡ (n -1)! is the Gamma function and v is a free variable. Replacing k = 2, to compute the contributions of the squares of X and Y results in:

where we have scaled the random variables by their standard deviation, σ n . Eq. 10, also known as the exponential distribution with scale parameter 1

We can transform Eq. 10 to obtain the distribution of the square-root of the addition of two random variables, as is done when taking the absolute value of a complex number pair (z) by applying the

The transformation of the free variable v → w(v) in the probability density function f (v) can be achieved by replacing the limits and applying the chain rule: We apply this relation to Eq. 10 to obtain the transformation v → w(v) by noting that w(v) = √ v, hence v = w 2 and dv dw = 2w. Eq. 10 then becomes:

Eq. 12 is also called the Gamma distribution. It has the interesting characteristic that its maximum-