Inference of Transmissivity in Crystalline Rock Using Flow Logs Under Steady‐State Pumping: Impact of Multiscale Heterogeneity

The significance of fracture roughness for hydraulic characterization of sparsely fractured (crystalline) rock using flow logs under steady‐state pumping conditions is investigated by numerical simulations in three‐dimensional fracture networks. A new solver for simulating flow in three‐dimensional discrete fracture networks is implemented using a hybrid finite volume and finite element method. The analysis is focused on different scales of heterogeneity in three‐dimensions: the network scale heterogeneity, the fracture‐to‐fracture scale heterogeneity, and the individual fracture scale heterogeneity (fracture roughness). The results show that the individual fracture scale heterogeneity due to fracture roughness is potentially an additional source of uncertainty for hydraulic tests in fractured rocks using flow logs. Specifically, it enhances the variation in measured flowrates and inferred transmissivity in addition to the network scale heterogeneity and the fracture‐to‐fracture scale heterogeneity. However, the individual fracture scale heterogeneity has relatively small impact on the inferred transmissivity distributions compared to the fracture‐to‐fracture scale heterogeneity. Moreover, fracture roughness and fracture‐to‐fracture heterogeneity have limited influence on the median values of the inferred transmissivity such that the median transmissivity is mostly reliable when inferred from flow logs. Comparison between the inferred and underlying input transmissivity distributions shows that interpreting hydraulic tests in crystalline rock using flow logs and the Thiem equation may underestimate the variation range of the underlying transmissivity. The results are helpful for understanding the uncertainty in hydraulic characterization of sparsely fractured rocks using flow logs that are relevant for various geo‐engineering and water resources applications.


Introduction
The hydraulic property of fracture transmissivity (i.e., the conductive property of fractures) directly controls groundwater flow and transport processes in a rock system since the fractures are the main pathways for mass and energy transport (Frampton & Cvetkovic, 2010;National Research Council, 1996;Neuman, 2005;Pickens et al., 1987), thus hydraulic characterization for the fracture transmissivity is important for many geo-engineering applications, including geological disposal of nuclear waste (e.g., Tsang et al., 2015), groundwater pollution control (e.g., Becker & Shapiro, 2003;Steimle, 2002), and geothermal extraction (e.g., Fox et al., 2016;Murphy et al., 1981). In particular, sparsely fractured crystalline rocks are considered as suitable geological formations for high-level radioactive nuclear waste disposal and hydraulic characterization of surrounding crystalline rocks is particularly important for the safety of long-term storage of nuclear waste (Bredehoeft & Maini, 1981;Tsang et al., 2015).
For crystalline rocks with low matrix porosity, rock fractures constitute the dominant flow conduits. The fracture networks control both overall permeability and individual flow/transport pathways (National Research Council, 1996). For hydraulic characterization of fracture transmissivity in crystalline rocks, pumping tests using high-resolution flow logs, especially the Posiva Flow Log (PFL) tool (Öhberg & Rouhiainen, 2000), have been widely used in site investigation for Finnish and Swedish radioactive waste disposal programs (e.g., Follin et al., 2007). However, hydraulic characterization of fracture transmissivity in crystalline rocks remains a challenge due to the strong heterogeneity of fracture networks with regard to both the fracture network structure and fracture transmissivity. Specifically, three scales of heterogeneity affect hydraulic properties of fracture networks: first, the inherent network scale (referred to as network heterogeneity) (de Dreuzy et al., 2012;Painter & Cvetkovic, 2001); second, the fracture-to-fracture scale (referred to as fracture-to-fracture heterogeneity) (Baghbanan & Jing, 2007;Hyman et al., 2016); and third, the individual fracture scale mainly caused by fracture surface roughness (also referred to as fracture internal heterogeneity) (e.g., Brown, 1987;Detwiler & Rajaram, 2007). The impact of these multiscale heterogeneities on flow and transport properties have been studied for the most part independently.
At the network and fracture-to-fracture scales, numerous studies applied discrete fracture network (DFN) models to investigate the hydraulic properties of crystalline rocks, both in 2-D (e.g., Baghbanan & Jing, 2007;Cacas et al., 1990;de Dreuzy et al., 2001;Long et al., 1982) and 3-D (e.g., de Dreuzy et al., 2012;Frampton & Cvetkovic, 2011;Hyman et al., 2015Hyman et al., , 2016Hyman et al., , 2019Maillot et al., 2016). In these studies, the transmissivity of each fracture is assumed constant in the network system, following assumed statistical distributions, such as power law (de Dreuzy et al., 2001;Hyman et al., 2016) or the lognormal distribution (Baghbanan & Jing, 2007). The randomly distributed transmissivity at the fracture-to-fracture scale could be uncorrelated, semicorrelated, or correlated to fracture size (Follin et al., 2007;Hyman et al., 2016;Neuman, 2008;Renshaw & Park, 1997). Using 3-D DFN modeling with consideration of the fracture-tofracture scale heterogeneity, Frampton and Cvetkovic (2010) discussed inference of field-scale transmissivity in crystalline rock based on PFL measurements.  inferred geometric and hydraulic parameters for a hydrogeological DFN model intended for site-scale groundwater flow and solute transport modeling in crystalline rock, based on detailed flow logging and down-hole digital imaging in cored boreholes. These studies have demonstrated that the high-resolution flow logs data can be very useful for inferring both geometrical and hydraulic properties of DFNs. However, the individual fracture scale heterogeneity was routinely neglected in these studies. Therefore, the impact of individual fracture scale heterogeneity on inference of fracture transmissivity using flow logs data remains unknown.
Because flow in sparsely fractured rock is controlled by conduits created by fractures, it is important to describe the aperture variations, from individual fractures to the network, in order to realistically simulate hydraulic tests. Aperture fields of natural rock fractures are spatially variable mainly due to fracture surface roughness (e.g., Brown, 1987;Detwiler & Rajaram, 2007). Extensive experimental and numerical studies have shown that the fracture roughness significantly affects the fracture transmissivity in single rock fractures (e.g., Brown, 1987;Brush & Thomson, 2003;Kang et al., 2016;Lee et al., 2014;Nicholl et al., 1999;Zimmerman & Bodvarsson, 1996;Zou et al., 2015Zou et al., , 2017a. In addition, a few experimental and numerical studies revealed that the local fracture roughness significantly affects channeling and mixing processes at rough-walled rock fracture intersections (Johnson et al., 2006;Li et al., 2020;Zou et al., 2017b). Moreover, a few numerical studies explored the impact of fracture roughness on the flow and transport processes in 3-D DFN systems (e.g., de Dreuzy et al., 2012;Frampton et al., 2018Frampton et al., , 2019Huang et al., 2019;Painter, 2006;Painter & Cvetkovic, 2001;Poteri & Laitinen, 1997). For instance, de Dreuzy et al. (2012) investigated the impact of fracture roughness on the equivalent permeability of a fractured medium using 3-D DFN simulations; it was found that the equivalent permeability of rough fracture networks is either larger (for dense fracture systems) or smaller (for sparse fracture systems) than the corresponding smooth fracture networks. Frampton et al. (2018Frampton et al. ( , 2019 investigated the impact of fracture roughness with connected and disconnected textures on advective transport in 3-D DFN; they found that the individual fracture scale heterogeneity may notably affect early mass arrival through a DFN. These studies generally revealed that the individual fracture scale heterogeneity is coupled with the network structures and the fracture-to-fracture heterogeneity, and hence will to some extent affect flow and transport properties of DFN systems. However, these studies have mostly focused on flow and/or transport processes under natural flow conditions, rather than pumping conditions for hydraulic characterization. The impact of individual fracture roughness on transmissivity inference in crystalline rocks using flow logs has not been comprehensively addressed. In particular, most interpretation models for hydraulic characterization of fractured rocks (e.g., Barker, 1988;Thiem, 1906) do not account for the impact of fracture roughness. To what extent the inferred transmissivity from pumping tests using flow logs (based on the assumption of neglecting fracture roughness) is representative of real systems remains an important open question.
The general objective of this study is to investigate the impact of multiscale heterogeneity including fracture roughness on transmissivity inference in crystalline rocks by flow logs under steady-state pumping. The specific objectives focus on answering the following questions: 1. What is the impact of multiscale heterogeneity on inferred transmissivity from flow logs under steady-state pumping in crystalline rocks? 2. To what extent is the transmissivity inferred from flow logs representative of the underlying transmissivity distributions, depending on multiscale heterogeneities?
To answer these questions, we simulate steady-state pumping tests in a 3-D DFN virtual setup, where the individual fracture aperture variability is modeled as a correlated lognormal random field in addition to the network and fracture-to-fracture heterogeneity. In practice, it is common to conduct multiple pumping tests at different spatial locations for a more comprehensive hydraulic characterization of a site. In this study, 20 cubic DFN models with the domain size of 100 m 3 based on the same site investigation data are adapted as being representative samples, to investigate the potential uncertainty among the sampling of pumping tests. The number of samples is selected according to sampling analysis. The result of sampling analysis is presented in the supporting information, showing that 20 samples can provide convergent statistical distribution of the measured flowrates.
The work is organized as follows. In the next section we first present the DFN virtual setup describing the multiscale heterogeneity that will be used for simulations and then summarize a hybrid method for solving water flow. In section 3, the results are presented and discussed. A Discussion section is then followed by a summary of main findings and conclusions.

DFNs Approach
Hydraulic tests using flow logs are modeled in 3-D DFN models where the fractures are represented as 2-D planes. The fractures are randomly generated by assigned distributions of location, shape, orientation, size, and density information that was obtained from site geological investigations. For the purpose of this study, the assigned distribution information is adopted from site investigations at Forsmark (Sweden), a potential site for Swedish radioactive nuclear waste disposal (e.g., Follin et al., 2007). The adopted distributions and corresponding parameters for generating the DFN model are summarized in Table 1. The parameter values Note. The values are adopted from Tables 11-20 (for the fracture domain FFM01 with the depth ranges from −200 to −400 m) in the SKB report R-07-48 (Follin et al., 2007). The expressions for the Fisher and lower law distributions can be found in Hyman et al. (2016).
are taken from those provided in Tables 11-20 in the SKB report R-07-48 (Follin et al., 2007). In addition, the shape of the fractures is assumed as square where the aspect ratio is 1. The location of the fracture (represented by coordinates of the fracture center) is assumed as following a uniform distribution.
By assigning fracture distribution and parameter information as listed in Table 1, the DFN model for this study is generated and meshed using the open source tool dfnWorks (Hyman et al., 2015), where a feature rejection algorithm is used for generating the 3-D DFN  and the LaGriT meshing toolbox (LaGriT, 2013) is used to generate conforming Delaunay triangular mesh for the DFN. The modules for fracture generation and triangular meshing are packed in the dfnWorks using Python scripts (Hyman et al., 2015). To simulate pumping tests, a 1-D pipe borehole with a constant radius is added in the center of the DFN model. The 1-D borehole is meshed by 1-D-line elements conforming to the DFN meshes at the intersecting points with corresponding fractures. Figure 1a shows one of the generated DFN models in a 100 × 100 × 100 m 3 domain with a uniform-sized mesh (around 0.2 m) for this study (see the supporting information). Note that the uniform-sized mesh is essential for considering the fracture roughness in the present study. Although the dfnWorks suite (Hyman et al., 2015) contains flow solvers using the parallelized subsurface flow and reactive transport code PFLOTRAN (Lichtner et al., 2015) and FEHM (Zyvoloski et al., 1988), we implemented a hybrid finite volume and finite element (HYFVE for short) method to solve flow in 3-D DFN models for convenience (Calgaro et al., 2008;Lal & Jabir, 2010).

Transmissivity Field With Fracture Roughness
The fracture-to-fracture heterogeneity for the transmissivity is often assumed as following the power law distribution (Follin et al., 2007) or the lognormal distribution (e.g., Baghbanan & Jing, 2007;de Dreuzy et al., 2001), supported or partially supported by experimental or field investigation data (Renshaw & Park, 1997). Many studies have shown that the transmissivity for each fracture may be (positively) correlated to its size (Follin et al., 2007;Renshaw & Park, 1997). Three types of correlation structures between fracture transmissivity and its size are often adopted in the literature, i.e., uncorrelated, semicorrelated, and perfectly correlated (Follin et al., 2007;Hyman et al., 2016). The impacts of such correlations on the flow and transport properties have been studied in the literature (e.g., Baghbanan & Jing, 2007;Follin et al., 2007;Hyman et al., 2016Hyman et al., , 2019Makedonska et al., 2016).
Without loss of generality, we consider that the fracture transmissivity is semicorrelated to fracture size by a power law relation, as the background fracture-to-fracture heterogeneity. Since we will focus on the impact of fracture internal aperture variability in the present study, the more idealized uncorrelated and perfectly correlated cases are not considered here. The fracture internal aperture variability is a source of uncertainty due to difficulties in field measurement of fracture roughness or the aperture fields. For theoretical and numerical studies, correlated multivariant Gaussian, lognormal, and fractal fields are often used to describe the fracture aperture variability (Frampton et al., 2019;Huang et al., 2019). In this study, we use the spatially correlated lognormal fields to represent the fracture internal aperture variability, which is in analogy to previous studies (e.g., Frampton et al., 2019) and is partially supported by laboratory measurement data for rock fracture apertures (e.g., Pyrak-Nolte et al., 1997). Therefore, the mean transmissivity of each discrete fracture is semicorrelated to its size and follows the power law, while the internal aperture field for each fracture is a

Water Resources Research
ZOU AND CVETKOVIC correlated lognormal field with assigned mean transmissivity (semicorrelated to its size), variance, and correlation length. The final transmissivity field for each fracture, T, combined with power law semicorrelated mean transmissivity and correlated lognormal fields for individual fractures can then be expressed as where L is the fracture size, a and b are the coefficients of the power law correlation with the fracture size, σ T is standard derivation of a lognormal distribution accounting for the uncertainty between fractures of the same size; N (0, 1) is a random variable sampled from a normal distribution with mean 0 and variance 1; σ is standard derivation of a lognormal field, and MG (0, 1, λ) are 2-D multi-Gaussian fields with the mean 0, variance 1, and correlation length λ. The first two terms represent the semicorrelated fracture transmissivity at the fracture-to-fracture scale and the third term represents the lognormal distributed fracture internal aperture variability.
The Fourier transform method is used to generate the 2-D correlated multi-Gaussian fields with zero mean and unit variance with the correlation length λ, given by (e.g., Hyman & Winter, 2014) where G(x,y) is a random uncorrelated Gaussian field; λ x and λ y are the correlation length in the x and y directions, respectively; FFT and IFFT represent the 2-D Fourier transform and inverse Fourier transform, respectively; and F(x,y) is the filter function, expressed as To demonstrate the impact of fracture internal aperture variability on hydraulic characterization of fractured rock using flow logs, a series of numerical simulations of pumping tests in the DFN models with different variances (i.e., σ 2 ¼ 0,1,and 3) of the lognormal fields for each fracture are conducted. Various correlation lengths (i.e., λ x ¼ λ y ¼ 0.75,1,1.25,and 1.5 m) of the lognormal fields were considered in preliminary testing but the impact was minor and the results are not presented here. The chosen ranges of these parameters are in analogy to those presented in the literature (e.g., Baghbanan & Jing, 2007;Frampton et al., 2019). σ 2 ¼ 0 represents the homogeneous case without fracture roughness, which is one of the two reference cases. The other reference case only considers the inherent network scale heterogeneity by assigning a constant transmissivity (equals to the mean fracture-to-fracture transmissivity in the homogeneous case when σ 2 ¼ 0) for all fractures, referred to as constant case. For all cases with consideration of fracture surface roughness, i.e., when σ 2 ¼ 1 and 3, the mean transmissivity for each fracture is assigned the same value to isolate the impact of network scale heterogeneity, where the semicorrelation parameters, i.e., a ¼ 1.6 × 10 −9 , b ¼ 0.8, and σ T ¼ 0.7, are inferred from the site investigation data (Follin et al., 2007).
For convenience, the 2-D multi-Gaussian fields are created on 2-D structured Cartesian grids with uniform size on a rectangular plane with dimensions larger than a considered fracture. The generated fields are then projected to the nodes of the triangular mesh for each fracture using barycentric interpolation. Figure 1 illustrates one of the 3-D DFN models with exemplified cases of fracture transmissivity variability.

Mathematical Model for Groundwater Flow in Rock Fractures
The mathematical model for groundwater flow in rock fractures is the Reynolds equation that simplifies the Navier-Stokes equation based on the lubrication approximation, i.e., assuming that the flow is laminar and fracture aperture is relatively small compared to its length (Zimmerman et al., 1991;Zimmerman & Bodvarsson, 1996). Let Ω be a fracture domain with a boundary where Γ D and Γ N denote portion of Γ prescribed with Dirichlet and Neumann boundary conditions, respectively. The mathematical model can then be written as where h D and q N are the prescribed boundary condition functions, and n is the outward normal vector of the boundary Γ N . For hydraulic analysis, the flowrate is of interest, which can be solved by the cubic law (Snow, 1965) where w is the fracture width, b is the fracture aperture, ρ is water density, g is gravity acceleration, and μ is water viscosity. Flow in the circular borehole can be calculated by (e.g., Bird et al., 1960) where r is the radius of the borehole. In practice, the boreholes can be up to 1,500 m long and the diameter ranges from 56 up to 200 mm. The typical measured borehole diameter is 76 mm (e.g., Sokolnicki & Heikkinen, 2008).
The flow equations in 3-D DFN models are solved using the implemented HYFVE method. This hybrid method is based on a node-centered control-volume scheme using a dual mesh system, where the finite volume mesh is built from the finite element mesh. It combines both advantages of finite volume and finite element methods, which preserves mass balance in local control volumes and directly provides nodal flowrate and velocity fields. The discretization scheme of this HYFVE method has been used in the field of computational fluid dynamics to solve the Navier-Stokes systems. Detailed information on the dual mesh system and the numerical scheme can be found in the literature (e.g., Calgaro et al., 2008;Lal & Jabir, 2010). The implemented HYFVE method has been verified against known analytical solutions for flow in single fractures.

Simulation and Results
For the simulation of pumping tests in the DFN model, the four vertical outer boundaries of the DFN model are assigned with constant head h ¼ 100 m, and the pumping water head at the top of the borehole located at the center of top boundary is h ¼ 90 m. The remaining top and bottom boundaries are set as no-flow, i.e., n · u ¼ 0. In the present study, packers are not considered, and the entire borehole is pumped. Exemplified results of simulated head distributions are presented in the supporting information.  (Figures 2a, 2c, 2e, and 2g) and the corresponding spatial flowrate distribution in the borehole (Figures 2b, 2d, 2f, 2h, and 2j) for one DFN sample. Only the fractures intersecting the borehole are shown here for a better visualization.

Distribution of Flowrates
In general, the distribution of the flowrates in both fractures and measured in the borehole are affected by conditions of hydraulic heterogeneity. For the first reference case with constant transmissivity, the flowrates exhibit relatively small variation in both discrete fractures (Figure 2a) and the borehole (Figure 2b). This variation in flowrate is caused purely by the network structure since the hydraulic property is set as constant in the whole DFN system for this case. The second reference case introduces fracture-to-fracture variation in hydraulic property at the network scale, where the transmissivity for each fracture is semicorrelated to its size. In this case (Figures 2c and 2d), the variation in flowrate becomes more significant compared to the constant case (Figures 2a and 2b), where the coefficient of variation (CV, the ratio between standard deviation and mean) is 3.5 times of that for the constant case. This result shows that fracture-to-fracture hydraulic heterogeneity significantly enhances the variation in measured flowrates.

Water Resources Research
For the cases with fracture roughness (Figures 2e-2h), the flowrate distribution in the fractures shows channeling flow patterns in each single fracture (Figures 2a and 2c). In particular, the channeling of flow becomes stronger when the variance of aperture fields increases from 1 (Figure 2e) to 3 (Figure 2h). This enhanced channeling flow in each single fractures is mainly caused by the fracture roughness, which is similar to results in previous studies for flow in single rough-walled fractures (e.g., Zou et al., 2017aZou et al., , 2017c or heterogeneous porous media (e.g., Gotovac et al., 2009). The fracture roughness further slightly enhances the variation in measured flowrates in the borehole. This result confirms that the fracture roughness will cause internal channeling of the flow and may further increase uncertainty for hydraulic tests in DFN systems in addition to the network structure and fracture-to-fracture hydraulic heterogeneity.
Empirical cumulative distribution functions (CDFs) and complementary CDFs (CCDFs) for all samples (the thin curves) with different cases of hydraulic heterogeneity are presented in Figure 3 to quantitatively illustrate the impact of different scales of hydraulic heterogeneity on measured flowrate distributions in numerical pumping tests. The ensemble (the thick solid curves) and average (the thick dash curves) CDFs and CCDFs for different scales of hydraulic heterogeneity are also plotted to illustrate the overall distribution features of all samples. In general, the slopes of the CDFs reduce and the tails of CCDFs become more significant with the increasing heterogeneity from the network scale to the combined internal fracture scale. The results of averaged CDFs and CCDFs confirm that the overall variability of measured flowrates by flow logs in the borehole becomes more significant with the semicorrelated network scale heterogeneity and the correlated log-normally distributed aperture fields. The increased variability of measured flowrates slightly increases with the increasing variance of aperture fields from σ 2 ¼ 1 (Figure 3c) to σ 2 ¼ 3 ( Figure 3d).

Inference of Fracture Transmissivity
One important task in hydraulic characterization is to infer the fracture transmissivity. In cases where pumping rates overwhelm natural flow, the fracture transmissivity can be inferred from the flowrate measured using flow logs in the borehole by applying the Thiem equation (Thiem, 1906), expressed by where r 1 is the radius of the borehole, r 2 is the radial distance to a boundary with constant head, Δh is a constant drawdown in the borehole in pumping test, and Q is the measured volumetric flowrate. In DFN modeling, the volumetric flowrate is calculated by multiplying the specific flowrate with the borehole perimeter. In practice, r 2 is an unknown variable and it is often assigned empirically, i.e., r 2 /r 1 ¼ 500 (e.g., Sokolnicki & Heikkinen, 2008). Nevertheless, it is an insensitive parameter since it appears in the logarithm term.
We calculate the fracture transmissivity using the Thiem Equation 9 based on the measured flowrates in numerical pumping tests in the DFN model and compare it with the assigned underlying transmissivity

Water Resources Research
ZOU AND CVETKOVIC distribution. By doing so, we could evaluate to what extent the inferred transmissivity is representative of the underlying transmissivity of the DFN system. Figure 4 shows comparison of the inferred transmissivity distributions with those of the input transmissivity for different scales of hydraulic heterogeneity, based on ensemble of the measured flowrates and input transmissivity for the 20 samples. For the constant case, although the constant input transmissivity is 1.28e-8 m 2 /s, the CDF and CCDF of the inferred transmissivity (the green curves) show that the inferred transmissivity is not constant but exhibits variation due to the network scale heterogeneity. By introducing the semicorrelated fracture-to-fracture variability of transmissivity to the constant case, i.e., when σ 2 ¼ 0, the CDF and CCDF of input transmissivity (the blue dash curves) represent the assigned semicorrelated fracture-to-fracture variability. The CCDFs of inferred and input transmissivity match well while the CDFs show discrepancy, indicating that the relatively large values of the inferred transmissivity match well with the input transmissivity and the relatively small values of the inferred transmissivity are smaller than that of the input transmissivity. With consideration of fracture roughness, i.e., when σ 2 ¼ 1 and 3, the CDFs and CCDFs of input transmissivity (the red and the black curves) show gradually increased ranges of variation with increasing σ 2 . In these cases, the CDFs of inferred and input transmissivity match well while the CCDFs show discrepancy, indicating that the relatively small values of the inferred transmissivity match well with the input transmissivity and the relatively large values of the inferred transmissivity are smaller than that of the input transmissivity.
To quantitatively analyze the variation in the samples and compare the inferred and underlying input transmissivity distributions for different scales of hydraulic heterogeneity, boxplot of the 20 samples for three different percentiles (tenth, fiftieth, and ninetieth) of the transmissivity are presented in Figure 5. These percentiles represent the relatively small (10%), the median (50%), and the relatively large (90%) values of the transmissivity. In general, these percentiles of the inferred transmissivity exhibit larger variability than those of the underlying input transmissivity.
For the distribution of the tenth percentile of the inferred transmissivity, the variation in the samples gradually increase in the cases with consideration of the fracture-to-fracture scale heterogeneity and the additional fracture internal heterogeneity. For the case that only considers the network scale heterogeneity (i.e., with constant transmissivity), the tenth percentile of the inferred transmissivity (the median is 4.5e-9 m 2 /s) is smaller than that of the underlying input transmissivity (where the input transmissivity is 1.28e-8 m 2 /s) (Figure 5a). For the cases with consideration of the additional scales of heterogeneity (σ 2 ¼ 0,1, and 3), the distributions of the tenth percentile of the inferred transmissivity are close to distributions of the underlying input transmissivity (Figure 5a), where the relative differences between the tenth percentile of the inferred and input transmissivity values are −31%, −8%, and 23% for the cases of σ 2 ¼ 0,1, and 3, respectively. This result indicates that considering fracture-to-fracture heterogeneity and fracture roughness in analysis of hydraulic tests using flow logs could improve the prediction of relatively small values of transmissivity.
The distribution of the median values of the inferred transmissivity distributions are generally close to the corresponding median values of the underlying input transmissivity for all cases (Figure 5b), where the relative differences between the median of the inferred and input transmissivity values are −29%, −25%, −27%, and −32% for the cases of constant, σ 2 ¼ 0,1, and 3, respectively. The distribution ranges of the median inferred transmissivity generally cover that of the underlying input transmissivity. This result indicates that hydraulic tests in sparsely fractured rocks using flow logs and the Thiem equation are robust and representative for inference of the median transmissivity. Meanwhile, the median values for all cases of hydraulic heterogeneity (even for the constant case) are nearly in the same order, i.e., around 1e-8 m 2 /s, which means that hydraulic heterogeneity both at the fracture-to-fracture scale and the individual fracture scale have limited influence on inference of the median transmissivity in hydraulic test using flow logs.
The distribution of the ninetieth percentile of the inferred transmissivity match well with that of the input transmissivity in the cases when fracture roughness is not considered, where the relative differences between the inferred and input transmissivity are 5% and 6% for the cases of constant and σ 2 ¼ 0, respectively. For rough fractures (σ 2 ¼ 1 and 3), the distribution of the ninetieth percentiles of the inferred transmissivity are generally below those of the underlying transmissivity (Figure 5c), where the relative differences between the inferred and input transmissivity are −79% and −206% for the cases of σ 2 ¼ 1 and 3, respectively. This means that the hydraulic tests of transmissivity using flow logs and the Thiem equation for sparsely fractured rocks underestimate the relatively larger values of the realistic transmissivity distributions. In addition, the distributions of the ninetieth percentiles of the inferred transmissivity for rough fractures are close to those for the case σ 2 ¼ 0. This indicates that the fracture roughness has limited impact on inference of the relatively larger values of transmissivity in hydraulic tests using flow logs compared to the impact of fracture-tofracture scale heterogeneity.
The overall variation ranges represented by the difference between the ninetieth and tenth percentiles for the inferred transmissivity are 8.9e-9, 3.7e-8, 4.2e-8, and 5.4e-8 m 2 /s for the four cases, respectively. The increasing variation ranges quantitatively illustrate the impact of different scales of heterogeneity on fracture transmissivity inference. It shows that considering fracture-to-fracture heterogeneity (σ 2 ¼ 0) significantly increases the overall variation range by 4.2 times compared to the constant case. Considering the additional internal heterogeneity of roughness further increases the overall variation range by 1.1 and 1.5 times compared to the case σ 2 ¼ 0, for the cases of σ 2 ¼ 1 and 3, respectively. For the underlying input transmissivity, the variation ranges are 0, 4.5e-8, 7.8e-8, and 1.7e-7 m 2 /s for the four cases, respectively. For the cases with consideration of the fracture-to-fracture heterogeneity and fracture roughness (σ 2 ¼ 1 and 3), the inferred transmissivity significantly underestimate the overall variation ranges, where the variation ranges of the inferred transmissivity are only 53% and 32% of the input transmissivity for the cases of σ 2 ¼ 1 and 3, respectively. Note that the absolute difference between the percentiles of the inferred and input transmissivity could be affected by using different parameter values in the Thiem equation, whereas the overall variation ranges presented here will not be affected. This means that hydraulic tests in sparsely fractured rocks using flow logs and the Thiem equation could underestimate variation ranges of the realistic transmissivity distributions.

Multiscale Heterogeneity
Our simulation results demonstrate impacts of different scales of hydraulic heterogeneity on transmissivity inference for crystalline rock using flow logs. Variations in the measured flowrates ( Figure 3) and the inferred transmissivity ( Figure 4) observed for the constant case indicate that the network structure is the instinct heterogeneity that affects the variability of the inferred transmissivity. The semicorrelated fracture-to-fracture heterogeneity enhances the variation in the measured flowrates and the inferred transmissivity. In addition to the fracture-to-fracture heterogeneity, the fracture roughness slightly magnifies the variation in the measured flowrates and the inferred transmissivity. As anticipated, transmissivity variations slightly increase with the variance of the fracture internal aperture variability. Our results show that the fracture-to-fracture heterogeneity is the primary scale of heterogeneity that affects transmissivity inference in addition to the instinct network scale heterogeneity. Compared to the fracture-to-fracture heterogeneity, the fracture internal heterogeneity of surface roughness has smaller impact on the transmissivity inference. This understanding is generally consistent with previous studies on fluid flow and solute transport in fracture networks under natural flow conditions (de Dreuzy et al., 2012;Frampton et al., 2018Frampton et al., , 2019Makedonska et al., 2016;Poteri & Laitinen, 1997).
At present, the impact of fracture internal heterogeneity of surface roughness is often ignored for transmissivity inference in practice, where only the mean transmissivities or the size correlation parameters are inferred from hydraulic test data through inverse modeling (e.g., Frampton & Cvetkovic, 2010). Our results confirm that ignoring the fracture surface roughness will have limited effect on transmissivity inference in practice. However, by ignoring fracture surface roughness, the inferred mean transmissivities or the size correlation parameters reflect the total impacts of multiscales hydraulic heterogeneity only in a statistical sense and they may not be consistent with the realistic (mechanical) DFN structures. Using such inferred hydraulic parameters and DFN structures to model solute transport or hydro-mechanical processes may introduce uncertainties which require further investigations (e.g., Frampton et al., 2019).
The median values of the inferred transmissivity using the Thiem equation (Figure 4) are mostly stable for all cases (e.g., the median values of the inferred transmissivity are of the same order, i.e., 1e-8 m 2 /s) and close to the median values of the underlying input transmissivity. This implies that both the fracture-to-fracture heterogeneity and the fracture roughness have limited influence on the median values of inferred transmissivity. Therefore, the median transmissivity inferred from pumping test using flow logs is representative for the underlying distribution and the flow log results can be effectively used for characterizing the median hydraulic properties for crystalline rocks. Moreover, the constant case by assigning the median of the inferenced transmissivity for a DFN model can be taken as the lower limit of variability for DFN modeling of flow and transport in fractured rocks. The upper limit of variability will be affected by the fracture-to-fracture variability and fracture roughness (e.g., de Dreuzy et al., 2012;Frampton et al., 2019;Hyman et al., 2016;Maillot et al., 2016).

Relating to Field Data
Crystalline rocks exhibit three dominant scales of heterogeneity that are addressed in this work: network scale, fracture-to-fracture scale, and the individual fracture scale. Therefore, the hydraulic test results obtained in practice reflect the uncertainty from the coupling of all three scales of heterogeneity. The impacts of different scales of heterogeneity on transmissivity inference shown in the present study are highly relevant for field pumping tests and can be used to interpret field data measured using high-resolution flow logs, such as PFL (Öhberg & Rouhiainen, 2000).
An example of a field pumping test data using the PFL conducted at the Laxemar subregion (it was one of the candidate repository sites for geological disposal of radioactive waste in Sweden) has been published in Frampton and Cvetkovic (2010), where their Figures 2 and 3 present the distributions of measured flowrates using PFL and the inferred transmissivity based on the Thiem equation, respectively. Our DFN modeling results (shown here in Figures 2 and 3) show similar distribution features to the field data presented in Frampton and Cvetkovic (2010) study (their Figures 2 and 3). Furthermore, the 3-D DFN model for simulating hydraulic tests with consideration of fracture internal scale heterogeneity could be used to infer the fracture surface roughness or aperture variability information from field pumping test data through inverse modeling.
We focused here on the impact of fracture roughness for a generic study and did not directly compare simulation results with field flow log data. In addition, we only considered the idealized steady-state pumping test conditions without considering any existing background natural flow. Inferring transmissivity and fracture surface roughness or aperture variability information from field pumping test data using flow logs by considering the more realistic boundary conditions under natural flow is an important topic for future study. Further investigations using the methodology proposed in this work may eventually improve the reliability of inferring different levels of heterogeneity in crystalline rock from a combination of hydraulic and structural field data (including flow logs) obtained, e.g., in the context of site characterizations .

Conclusions
In this work, we investigated the impact of multiscale heterogeneity on hydraulic characterization of crystalline rocks using 3-D DFN modeling that utilizes a HYFVE flow solver. The most important conclusions from this study are summarized as follows: 1. Essentially three scales of heterogeneity control flow in crystalline rocks in a 100 m block: network scale, fracture-to-fracture scale, and fracture internal scale (fracture roughness). The network scale is the instinct heterogeneity for the observed variations of flowrate and the inference of transmissivity from flow logs under steady-state pumping. These variations will be significantly enhanced by the fractureto-fracture scale heterogeneity, where the overall variation range (represented by the differences between the ninetieth and tenth percentiles) of the inferred transmissivity increases 4.2 times compared to the case only considering the network scale heterogeneity. 2. Fracture roughness further enhances variation in the measured flowrate and the inferred transmissivity.
The enhanced variation increases with the increasing fracture roughness variance. However, compared to the fracture-to-fracture heterogeneity, the fracture roughness has smaller impact on the inferred transmissivity distributions, where the overall variation range of the inferred transmissivity increases 1.1 and 1.5 times for the cases with fracture roughness variance σ 2 ¼ 1 and 3, respectively. 3. Both the fracture-to-fracture heterogeneity and internal heterogeneity have limited impact on the inferred median transmissivity under steady-state pumping, where the median transmissivity obtained from the Thiem equation applied on observed flowrates is close to the median of the underlying transmissivity distribution, where the maximum relative difference is less than 32%. Therefore, the inferred median transmissivity is a robust measure of hydraulic properties that can be used to determine the lower limit of hydraulic variations in a fractured rock system. 4. Beside of the median values, the inferred transmissivity from the Thiem equation using flow logs (under steady-state pumping) generally underestimate the overall variation range of the underlying realistic transmissivity. Specifically, the variation ranges of the inferred transmissivity are only 53% and 32% of the input transmissivity for the cases with fracture roughness variance σ 2 ¼ 1 and 3, respectively.
Directly measuring aperture fields for large-scale natural rock fractures still poses considerable challenges (e.g., Detwiler et al., 1999;Pyrak-Nolte et al., 1997). Fracture apertures have been assumed as generic random fields, correlated or uncorrelated (e.g., de Dreuzy et al., 2012;Frampton et al., 2019;Huang et al., 2019;Makedonska et al., 2016). The impact of more realistic fracture roughness and aperture structures supported by high-resolution laboratory or field data, on transmissivity inference for crystalline rocks using flow logs, need to be further studied. We considered here the idealized steady-state pumping test using flow logs with fluid flow as the single physical process. The impact of fracture roughness on transient hydraulic tests with hydro-mechanical coupling, and comparison with field data in settings that combine natural flow and pumping, are topics to be investigated in the future. In addition, the potential underestimation of variation ranges of the inferred transmissivity in hydraulic tests using flow logs and the Thiem equation demonstrated in the present study and its effects on modeling of mass and energy transport need to be further investigated with consideration of the more realistic conditions in the field.