Conditioning of stochastic 3-D fracture networks to hydrological and geophysical data Caroline Dorna, Niklas Lindea, Tanguy Le Borgneb, Olivier Bourb, 7 Jean-Raynald de Dreuzyb 1 0 2 aUniversity of Lausanne, Faculty of Geoscience and the Environment, Applied and n Environmental Geophysics Group, Switzerland a bUniversit´e Rennes 1-CNRS, OSUR G´eoscience Rennes, Campus Beaulieu, 35042 J Rennes, France 7 ] h p - Abstract o e g Thegeometryandconnectivityoffracturesexertastronginfluenceonthe . s c flow and transport properties of fracture networks. We present a novel ap- i s y proach to stochastically generate three-dimensional discrete networks of con- h p nected fractures that are conditioned to hydrological and geophysical data. [ 1 A hierarchical rejection sampling algorithm is used to draw realizations from v 1 theposteriorprobabilitydensityfunctionatdifferentconditioninglevels. The 8 8 method is applied to a well-studied granitic formation using data acquired 1 0 within two boreholes located 6 m apart. The prior models include 27 frac- . 1 0 tures with their geometry (position and orientation) bounded by informa- 7 1 tion derived from single-hole ground-penetrating radar (GPR) data acquired : v i during saline tracer tests and optical televiewer logs. Eleven cross-hole hy- X r draulic connections between fractures in neighboring boreholes and the order a in which the tracer arrives at different fractures are used for conditioning. Furthermore, the networks are conditioned to the observed relative hydraulic importance of the different hydraulic connections by numerically simulating Preprint submitted to Advances in Water Resources January 10, 2017 the flow response. Among the conditioning data considered, constraints on the relative flow contributions were the most effective in determining the variability among the network realizations. Nevertheless, we find that the posterior model space is strongly determined by the imposed prior bounds. Strong prior bounds were derived from GPR measurements and helped to make the approach computationally feasible. We analyze a set of 230 pos- terior realizations that reproduce all data given their uncertainties assuming the same uniform transmissivity in all fractures. The posterior models pro- vide valuable statistics on length scales and density of connected fractures, as well as their connectivity. In an additional analysis, effective transmissivity estimates of the posterior realizations indicate a strong influence of the DFN structure, in that it induces large variations of equivalent transmissivities between realizations. The transmissivity estimates agree well with previous estimates at the site based on pumping, flowmeter and temperature data. Keywords: discrete fracture network, conditioned fracture network, data integration, ground-penetrating radar, probabilistic inversion 1. Introduction The characterization of flow, storage, and transport properties in frac- tured rock formations is challenging (e.g., Bonnet et al., 2001; Berkowitz, 2002; Neuman, 2005). The challenge resides mainly in the extreme hetero- geneity of fractured formations, in which hydraulic conductivity varies over many orders of magnitudes within very short distances. Furthermore, flow and transport in fractured systems are highly organized, with flow channel- ing within fracture planes and preferential pathways within the connected 2 fracture network (e.g., Smith and Schwartz, 1984; Neretnieks, 1993). Di- rect measurements of hydraulic or geometrical properties of all individual fractures are impossible, and the data coverage in field investigations is typ- ically very low compared with the underlying heterogeneity. One solution to this problem is to derive effective bulk hydraulic properties (Mar´echal et al., 2004), but this approach has been widely criticized, as it implies the exis- tence of a homogenization scale that is smaller than the scale of investigation (Berkowitz, 2002; Neuman, 2005). Ongoing research focuses on the hydraulic and geometrical characterization of fracture networks (e.g., network perme- ability, fracture density, fracture length distribution), of individual fractures (e.g., position, location, orientation, transmissivity) and their connectivity (e.g., Berkowitz, 2002). Hydrological investigations of fractured rock aquifers often include time- consuming hydraulic, flowmeter and tracer tests (e.g., Harvey and Gorelick, 2000; Day-Lewis et al., 2006; Le Borgne et al., 2007). These tests allow in- vestigating hydraulic properties and dispersive phenomena in transmissive fracture networks that intersect the boreholes. Flowmeter and tracer tests provide effective parameter estimates for the investigated hydraulic connec- tions, which might comprise many fractures with complex connectivity pat- terns (e.g., Frampton and Cvetkovic, 2010). Inferring hydraulic characteris- tics of the individual fractures is most often impossible, unless if conditions arefavorable(Novakowskietal.,1985). Evenextensivepackertestsfailtore- solve geometrical and hydraulic properties of individual fractures (Hao et al., 2007). Previous geophysical field investigations in fractured rocks have demon- 3 strated the value added of single-hole ground-penetrating radar (GPR) re- flection data in providing constraints on the geometry of individual fractures (e.g., Olsson et al., 1992; Dorn et al., 2012a). Properties, such as their loca- tion and orientation relative to the borehole, can be constrained from GPR images (e.g., Olsson et al., 1992; Grasmueck, 1996; Tsoflias et al., 2004). Compared with other geophysical methods that are used in fractured rock investigations, GPR is arguably the most favorable tool to image individual fractures away from boreholes. Seismic methods, for instance, typically em- ploy sources with a resulting larger wavelength and the acoustic impedance contrast is less pronounced than the electromagnetic impedance in most frac- tured rock environments (e.g., Mair and Green, 1981; Khalil et al., 1993). The resolution of electrical resistivity tomography is insufficient to identify individual fractures (Day-Lewis et al., 2005; Robinson et al., 2013; Slater et al., 1997), but this is also the case for many other types of geophysical data when interpreted by smoothness-constrained inversion (such as seismic or GPR attenuation or travel time tomography; e.g., Day-Lewis et al. (2003); Ramirez and Lytle (1986); Daily and Ramirez (1989)). Surface-deployed GPR can be used to image shallow dipping fractures (e.g., Grasmueck, 1996), whereas vertical borehole acquisitions can be used to image steeply dipping fractures (e.g., Olsson et al., 1992). Borehole GPR reflection data can either be acquired in single-hole or cross-hole mode (Dorn et al., 2012a). Common borehole GPR systems create omni-directional radi- ation patterns (symmetry around the borehole axis), which precludes the de- terminationofazimuthsofplanereflectorswhenusingdatafromoneborehole alone. Directional antennas that were originally developed for the character- 4 ization of prospective nuclear waste repositories, have only recently become more widely developed, but are still not commonly used (Slob et al., 2010). The deployment of GPR monitoring during saline tracer tests is useful to identify those fractures that are transmissive and connected to the injection point. Tracer movements can thus be imaged in individual fractures away from the borehole locations (e.g., Talley et al., 2005; Dorn et al., 2012b). In the following, we refer to this data type as time-lapse GPR data. The accessibility of fractured rock formations are typically restricted to the boreholes. This restriction has important consequences, in that the ori- entation of fractures can be well resolved at the borehole locations (using optical or sonic televiewer logs), but only within significant uncertainty away from the borehole (e.g., using GPR reflection images) (Olsson et al., 1992; Dorn et al., 2012a). Furthermore, hydrological and geophysical data gener- ally constrain mainly hydraulic or geometrical properties, respectively. The integration of different data types can help to overcome some of the limita- tions commonly encountered when characterizing fracture networks. As both hydrological and geophysical data and their interpretations contain consid- erable uncertainties, they should be carefully taken into account within a formal data integration procedure (e.g., Le Goc et al., 2010). The use of geophysical data in building aquifer models has been investi- gated thoroughly in recent years (e.g., Hubbard et al., 1997; Gallardo and Meju, 2007; Linde et al., 2006; Day-Lewis et al., 2000; Chen et al., 2006). Many successful field demonstrations of deterministic integration methods exist that determine hydraulic and structural properties of unconsolidated materials (e.g., Looms et al., 2008; Doetsch et al., 2010), but fewer exam- 5 ples consider fractured rock systems (e.g., Day-Lewis et al., 2003; Linde and Pedersen, 2004). Deterministic inversion algorithms seek one model that fits the data, while deviating the least from a reference model or preferred model morphology. Their applicability may be limited for highly non-linear problems as the typically used gradient-based optimization techniques may get trapped in local minima and they provide limited insights about pa- rameter uncertainty (e.g., Tarantola, 2005). Probabilistic inversion methods are suitable when it is expected that very different models (e.g., in terms of fracture geometries) are consistent with the available prior information and data, as they may explore the posterior probability density functions (pdfs) (e.g., Mosegaard and Tarantola, 1995). Probabilistic sampling-based approaches for high-dimensional problems (e.g., hundreds of parameters) are computationally expensive. This means that it is often necessary to decrease parameter dimensionality and their ranges as much as possible beforehand. Numerous global search methods exist, for example Markov chain Monte Carlo methods (MCMC), simulated annealing or genetic algorithms. Day- Lewis et al. (2000) conditioned the 3-D geometry of fracture zones to hy- draulic data using a simulated annealing approach. MCMC has the advan- tage that parameter uncertainty is formally assessed (e.g., Hastings, 1970; Mosegaard and Tarantola, 1995; Laloy and Vrugt, 2012). Chen et al. (2006) used MCMC sampling to estimate probabilities of zones having high hy- draulic conductivities. Markov chain Monte Carlo methods generate param- eter samples that converge to a stationary distribution that coincides with the posterior pdf under certain conditions (e.g., Hassan et al., 2009; Chen et al., 2006). Such methods can be inefficient for discrete fracture networks 6 (DFN) for which small changes in fracture geometries might lead to very dif- ferent model responses. In this case it may be efficient to use the brute force rejection sampling method, which is the only exact and also the simplest global search method (Mosegaard and Tarantola, 1995). To investigate the 3-D heterogeneity of fracture networks, the stochastic generation of multiple DFN models (e.g., Darcel et al., 2003) that are con- ditioned to available data is appealing. The conditioning of DFN models to borehole observations is common, but its utility has been questioned given the difficulty in correlating fracture aperture at the borehole locations to hydraulic apertures (Renshaw, 1995) and the limited coverage of borehole data with respect to the investigated rock volume (Neuman, 2005). DFN approaches rely on statistical descriptions of fracture distributions, but reli- able statistics are often very difficult to obtain. To partly circumvent these problems, we propose herein to condition connected 3-D DFN models to hy- draulic, tracer and GPR reflection data that are sensitive to fracture geome- try and hydrological state variables. Also, the flow responses of the realized DFNs are numerically simulated. For the conditioning of the DFNs we use a hierarchical rejection sampling method. We apply our scheme to data from the Ploemeur site in France, where relatively few well-connected fractures appear to control the flow (Dorn et al., 2012b). More specifically, we derive 3-D models of connected DFNs by combining (1) logging data together with (2) a few fundamental geometrical properties of connected fracture networks, (3) images of fractures obtained from GPR data acquired under natural con- ditions and during tracer test, (4) topological constraints that describe the order in which the tracer arrives in different fractures as inferred from time- 7 lapse GPR data, and (5) hydrological information (e.g., tracer arrival times, massflux curves) inferred from tracer tests. The set of conditioned fracture models can then be used to evaluate fracture network statistics. 2. Methodology Within this study, we generate and condition connected DFN models un- der the assumption of an impermeable rock matrix. Fractures are described as homogenous thin circular discs distributed in 3-D space, where each frac- ture is parameterized with a midpoint (x,y,z), a dip (ϕ), an azimuth (θ) and a radius (R). The same uniform transmissivity is assigned to every frac- ture. The geometry of all individual fractures in a given model is defined by a vector m. A series of observed hydrological and geophysical data d obs stems from the noise-contaminated response of a true and unknown system. The goal of the probabilistic data integration or joint inversion procedure is to draw samples from a pdf p(m) describing prior constraints on the model space that also agree with d within the associated data and modeling er- obs rors. The agreement between calculated data d and observed data d cal obs is evaluated by a likelihood function L(m|d ). Bayes theorem is used to obs calculate the posterior pdf p(m|d ) given the prior p(m) and the likelihood obs function L(m|d ), where L(m|d ) equals p(d |m), obs obs obs p(m|d ) = const p(m)L(m|d ) (1) obs obs where const is a normalization constant (e.g., Mosegaard and Tarantola, 1995). Many different methods have been developed to derive p(m|d ). obs Rejection sampling is the only exact Monte Carlo sampling technique and 8 every accepted model is a random and independent draw from p(m|d ). obs This simple method draws proposals directly from the prior distribution and their acceptance is decided by the value of the likelihood function alone, that is, there is no comparison with the likelihoods of previously sampled models, as in MCMC. A suitable metric is needed to compare simulated and observed data for a proposed model. Depending on the data type, we use the infinity norm l or ∞ the L2-norm l norm. The l norm is applicable if the error bounds are strict 2 ∞ (Marjoram et al., 2003). In this case, a model is accepted if all the residuals are within the error bounds of the observed data. The likelihood function is thus either zero or one. It may have complex patterns for fracture networks whose connectivity may change significantly with only small perturbations in fracture geometry. Different simulations are used to evaluate the model responses of each data type and the associated computational costs are typically very different (e.g., flow simulations are computationally more expensive than evaluating hydraulic connections). To keep the computational costs as low as possible, we apply a hierarchical formulation of Bayes theorem (e.g., Glaser et al., 2004), in which the samples from the posterior pdf at one conditioning level are subsequently evaluated at the next level using other data types. Time- consuming simulations of proposed models are thus only performed on those that agree with the data types for which quick simulations are available. Rejection sampling can be very slow when dealing with large data sets with small errors as the acceptance rate may be prohibitively low. This algorithm applies therefore primarily to instances when there are compara- 9 tively few data and when data or modeling errors are rather large. We argue that this is often the case in fractured rock investigations. Furthermore, comparatively narrow bounds on the prior distribution are crucial to make the rejection sampling approach feasible. For example, it is most useful to know how many fractures are part of the connected DFN and their approx- imate positions. To make this approach feasible, it is therefore necessary to use site-specific data to derive a bounded p(m) (these data cannot be used further as conditioning data in the likelihood function). Outmost care is needed to assure that these bounded priors reflect realistic distributions, thereby, avoiding over-confident predictions or incompatibilities between the prior and the likelihood terms. In this work, the properties of individual fractures are largely constrained by the prior bounds and subsequent condi- tioning is mainly used to assure that the proposed models are in agreement with data that characterize the behavior of the whole fracture network (e.g., hydraulic cross-hole connections) that are difficult to include in p(m). To demonstrate the methodology, we consider a set of classical hydrolog- ical data: optical televiewer logs, flowmeter and tracer test data. Further- more, we use single-hole GPR data that provide geometrical constraints on individual fractures that are part of the connected fracture network (Dorn et al., 2011, 2012b). Indeed, the method capitalizes through its model pa- rameterization on a detailed description of the fracture geometry derived from GPR imagery and televiewer data. The different data types and their uncertainties are described in the following. (1) Optical televiewer logs constrain the intersection depth, dip and az- imuth of fractures that intersect the boreholes. The uncertainties of these 10