Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

The effect of spatial aperture variations on the thermal performance of discretely fractured geothermal reservoirs

The effect of spatial aperture variations on the thermal performance of discretely fractured... cornell.edu Cornell Energy Institute, The effect of spatial aperture variations on the thermal performance of discretely frac- Cornell University, tured geothermal reservoirs was investigated using finite element method solutions to Ithaca 14850, USA the convective heat transfer in the fracture coupled with a boundary integral descrip- Full list of author information is available at the end of the tion of conductive heat transfer in the rock. The dipolar flow was generated by a source article of fluid volume at an injection well and a sink at a production well within a circular fracture. The statistics of the thermal performance of an ensemble of fracture realiza- tions was evaluated. Fractures with self-affine aperture fields where long range correla- tions dominant over short range correlations were generated. The results showed that spatial aperture variations most frequently lead to diminished thermal performance by creating flow channeling that reduces the heat transfer area. Enhanced thermal performance occurred in the less common cases when the aperture was small in the region between the wellbores, causing fluid flow to sweep out greater areas of the fracture and extract heat from a larger area. The standard deviation of the apertures had the largest influence on the thermal performance, while the spatial correlation of the aperture played a secondary role. Larger values of the standard deviation led to more adverse thermal performance. For the range of standard deviations investigated, the fraction of fractures exhibiting enhanced thermal performance compared to the base case of no aperture variations ranged from 34 to 49 %. The degradation of thermal performance due to aperture variations was largest when the well bore spacing was a larger fraction of the fracture diameter. Reservoirs consisting of two non-intersecting fractures connected to the same injection and production wells were also modeled. The uneven split of the flow between the reservoirs in this case caused a further deterioration of thermal performance compared with the single fracture reservoirs. However, flow control was able to overcome nearly all of the additional loss of thermal performance for the multifracture reservoir. Keywords: Enhanced geothermal systems, Aperture variations, Roughness, Discrete fractures, Heat transport in fractures , Coupled convective conductive heat transport, Thermal hydraulic modeling, Geothermal reservoir performance modeling Background Variations in fluid velocity on various flow paths in discretely fractured geothermal res - ervoirs like those found in Enhanced/Engineered Geothermal Systems (EGS) have a sig- nificant effect on the thermal performance of a reservoir. One source of variations in © 2015 Fox et al. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Fox et al. Geotherm Energy (2015) 3:21 Page 2 of 29 fluid flow is from spatial aperture variations in the fractures. EGS reservoirs are devel - oped in hot rock with low permeability and porosity to extract a portion of the stored thermal energy contained in the rock mass. A reservoir is stimulated by pressuriza- tion to open existing sealed fractures (hydro-shearing) or induce new ones by hydraulic fracturing where water may be pumped into and out of the fractured reservoir using a connected set of injection and production wells to “farm” or extract the stored ther- mal energy. While convective heat transport is the dominant mode of transport within the fractures, it is inherently coupled to conductive heat transport from the surrounding −6 2 bulk rock. Because of the low thermal diffusivity of the host rock (of order 10 m /s), only a small region surrounding the fracture network is effective for heat extraction. Given the finite areal extent of the fracture network, extended use will locally cool down the reservoir. The mass loading per unit area usually expressed as the total flow rate per unit fracture surface area is a critical parameter that often scales directly with the rate of thermal drawdown. Drawdown is a major concern for geothermal developers who must forecast reservoir performance and consider its effect in their techno-economic analysis. This study seeks to determine the effect of fracture aperture variations on the thermal performance of discretely fractured geothermal reservoirs Because fractures do not have smooth surfaces, spatial aperture variations can lead to maldistributed flow within the fracture network adding a degree of thermal performance uncertainty. Fracture core samples have revealed fracture apertures to have a self-aff - ine spatial correlation governed by a roughness or Hurst coefficient (Brown et al. 1986; Glover et al. 1998; Schmittbuhl et al. 1993, 1995; Plouraboué et al. 1995; Boffa et al. 1999; Ponson et  al. 2007). A self-affine spatial correlation is one characterized by long range correlations. The aperture undulations in the core samples are found on all length scales but the aperture variance is dominated by the largest scale. Often self-affine fracture surface variations are referred to as fracture roughness despite the fact that roughness implies variations on small rather than large length scales. In this paper, we have avoided using “rough and roughness” as descriptors as much as possible and instead use terms like aperture or surface variations. To determine the effect fracture surface variations have on fluid flow and on a reser - voir’s thermal performance, a dipole flow in circular fractures between a single injector and producer was selected. Dipole flow is not uniform and is characterized by stream - lines where the fluid velocity is faster near the injection and production wells. For sim - plicity, effects like temperature dependent properties and buoyancy were not considered. We are interested in modeling reservoirs with different degrees of fracture surface varia - tions, roughness/Hurst coefficients, wellbore spacings, and multiple fractures. Our focus is to determine how fracture aperture variations affect reservoir parameters like the pro - duction temperature and the likelihood a reservoir will exhibit thermal enhancement or decline through statistical parameters involving an ensemble set of fracture realizations. Through these statistical values, we hope to 1. infer how the extent of aperture variations aid or hinder thermal performance, 2. characterize the degree of reservoir performance uncertainty due to aperture variations, and 3. determine under what scenarios aperture variations need to be modeled. Furthermore, modeling spatial aperture variations pro- vides a method to systematically generate variations in flow fluid rates and their pathway through the fracture system. Fox et al. Geotherm Energy (2015) 3:21 Page 3 of 29 While little work has been done on the coupled convective-conductive heat trans- port of the fracture and surrounding rock in a reservoir with spatial aperture variations, researchers have modeled fluid flow and other transport properties using a self-affine aperture field. They have observed that self-affine aperture fields lead to flow chan - neling where most of the fluid will travel through several preferential channels (Brown 1987; Tsang and Tsang 1989; Glover et al. 1998; Drazer and Koplik 2000; Méheust and Schmittbuhl 2001, 2003; Auradou et  al. 2006). Méheust and Schmittbuhl (2001) deter- mined that three parameters are required to model the flow, the mean aperture, the root mean square of the aperture, and the Hurst coefficient. They further showed how the mean flow rate was either enhanced or inhibited, depending on the orientation of the applied pressure drop. Transport properties such as dispersion have also been studied in self-affine fractures (Plouraboué et al. 1998; Drazer and Koplik 2002; McDermott et al. 2009; Auradou et al. 2010). Numerical studies exhibit non-Gaussian transit time distri- bution with long tails (Drazer and Koplik 2002) and experiments with Newtonian and shear thinning fluids showed that geometric/hydraulic dispersion is dominant at low Péclet numbers with the dispersion coefficient on the order of the Péclet number (Aura - dou et al. 2010). Other researchers have used flow channeling to reduce numerical computation by reducing flow in 2D fractures to an equivalent 1D flow (Cacas et al. 1990; Nordqvist et al. 1992; Bruel et  al. 1994) or to reduce 3D flow into an equivalent 2D system (Bruderer- Weng et al. 2004). The above mentioned studies have focused on the effects fracture sur - face variations have on hydraulic performance and particle/solute transport and do not cover the consequence on convective heat transport. However, these studies mention preferential flow paths which suggest that surface variations may have a negative impact on thermal exchange. Kolditz (1995) suggested that local aperture variations may play a strong role in a res- ervoir’s thermal performance. Jupe et al. (1995) attempted to model spatial aperture var- iations by modeling separate flow paths. They used the multifracture parallel rectangular reservoir model described by Nicol and Robinson (1990) and treated the different flow paths as separate fractures. Their work focused on determining the wellbore spacing needed to obtain an economic EGS reservoir. By incorporating multiple flow paths, the wellbore spacing needed to obtain their desired performance criteria went from 700 m (no channeling) to 1080 m. The uncertainty of reservoir parameters like heat capacity and permeability have been modeled for a heterogeneous porous media geothermal system (Watanabe et al. 2010). The spatial correlation of the reservoir parameters were modeled using a spherical vari - ogram. Using a hundred realizations of the reservoir, Watanabe et al. (2010) determined that permeability was the most sensitivity parameter, with the production temperature having an uncertainty range of 40 K for the 15 years. Since the aperture can be viewed as an analog to permeability, we would expect spatial aperture variations to substantially control a reservoir thermal performance. The thermal effect of surface variations has recently been investigated in a rectangular fracture where the base flow is uniform (Neuville et al. 2010a , b). However, Neuville et al. (2010a, b) used a simplified heat transfer model that only considers the heat transfer in the fracture and not in the surrounding rock; the surrounding rock was assumed to be Fox et al. Geotherm Energy (2015) 3:21 Page 4 of 29 permanently at the original reservoir temperature. As a consequence, their results can be interpreted as the early thermal behavior of the reservoir and not the behavior after extended reservoir use. The assumption of constant rock temperature is valid for the time scale of w /α, where w is the fracture aperture and α is the thermal diffusivity of the −6 2 rock. For an aperture value from 1 to 5 mm and a thermal diffusivity of 1 × 10  m /s, the time scale will be on the order of 1–10  s. Neuville et  al. (2010b) work focused on developing relationships between a reservoir’s hydraulic performance and its thermal performance. They concluded heat exchange in self-affine fractures is less efficient com - pared to a constant aperture case. As a case study, Neuville et  al. (2010a) modeled the Soultz-sous-Forêts EGS reservoir as a single fracture and concluded that merely know- ing the mean aperture is not sufficient to determine the thermal heat exchange and that fracture surface variations always led to inhibited thermal exchange. Although Neuville et al. (2010a, b) have modeled variations in convective heat transport in the presence of aperture undulations, there is still a need to model the coupled convective heat trans- port in the fracture with that of the conductive transport in the surrounding rock. Our approach fully couples the heat transfer in the fracture to the surrounding bulk rock by using a finite element method in the fracture and boundary integral equations to describe the temperature in the bulk rock. Additionally, our study aims to determine how much deviation in production temperature in the absence of aperture undulations is anticipated after prolonged reservoir utilization. The use of statistical values of an ensemble of realizations including the mean and standard deviation of production tem- peratures will offer a better and more practical understanding of the effects fracture sur - face variations have on discretely fractured geothermal reservoirs. Methods The following section outlines our numerical approach in modeling self-affine aperture fields and the thermal performance of the resulting reservoirs. We first introduce the mathematical theory of self-affine aperture fields and describe the algorithm used to generate the fracture surfaces. Then, we detail the method used to solve for the fluid flow in the fracture and finally, we describe our numerical approach in solving for the heat transport in the reservoir and the transient temperature field. Self‑affine aperture field A spatially varying aperture field fluctuates with a value of h(x,  y) from a mean value of w to give w(x, y) = �w� + h(x, y). (1) The spatial mean aperture w is given by w(x, y) dA �w� ≡ , (2) dA where A is the fracture surface area and the variance or mean square of the aperture field is given by h (x, y) dA 2 2 A σ ≡ h = . dA A (3) Fox et al. Geotherm Energy (2015) 3:21 Page 5 of 29 The spatial correlation of an aperture field can be described by the variogram function γ, which is twice the difference of the aperture variance and the two-position autocor - relation function: ′ ′ 2 ′ γ(r ) ≡ h(r) − h(r + r ) = 2 h − h(r)h(r + r ) , (4) where r is the displacement from point r and r is a point on the surface of the fracture. The autocorrelation function is a measure of the correlation of the aperture at point r with that at point r + r . If the variogram function follows ′ ′ 2ζ γ(r ) = C|r | , (5) ′ 2ζ ′ it has a scaling behavior of γ(r ) =  γ(r ), where  is a scaling parameter and the field is said to have a self-affine geometry governed by the Hurst or roughness coefficient, ζ (Feder 1988; Falconer 2003). For this study, we have chosen to model fracture aperture variations as having a self-affine geometry. By measuring a surface profile, many mate - rials have been shown to have a self-affine geometry (Bouchaud 1997). From fracture core samples, the Hurst coefficient has been determined to be from 0.7 to 0.8 (Brown et al. 1986; Glover et al. 1998) and, specifically, equal to 0.8 for granite (Schmittbuhl et al. 1993, 1995; Plouraboué et  al. 1995), the typical rock type of deep geothermal systems. Sandstone, on the other hand, has been shown to have a lower value of ζ at around 0.5 (Boffa et  al. 1999; Ponson et  al. 2007). Although it is difficult to determine whether a fracture will exhibit a self-affine behavior at the length scale of geothermal fractures, it has been observed up to the meter scale (Brown and Scholz 1985). With a self-affine behavior, larger scale variations in h (small frequencies, long wave - lengths) are more important than small scale (large frequencies, short wavelengths) vari- ations. The large scale variations become more strongly dominant when ζ is increased. However, the large scale variations are the primary determinant of the aperture varia- tions for any positive ζ. A self-affine aperture field is analogous to fractional Brownian motion, if one views the spatial displacement |r | as a time-like variable and the change of aperture with displacement as equivalent to the spatial step of the Brownian motion. ζ< 1/2 With this analogy, corresponds to a subdiffusive or antipersistent process, ζ = 1/2 corresponds to a normal diffusive process, and 1/2 <ζ < 1 corresponds to a super diffusive process with ζ = 1 representing ballistic motion (Mandelbrot and Ness 1968). Interested readers are advised to consult the following references for more details on the theory of self-affine behavior and other applications, (Feder 1988; Peitgen and Saupe 1988; Falconer 2003). The algorithm used to generate spatially varying aperture fields, which was adapted from Méheust and Schmittbuhl (2001), uses random seeding and Fast Fourier Trans- forms to obtain isotropic self-affine fracture aperture fields on rectangular fractures. In order to generate circular fractures with a self-affine aperture field, a square fracture of length 2R was first generated and only the area of an inscribed circle was used for the analysis. Any shaped fracture could have been generated by first generating a large enough rectangular fracture to inscribe the desired shape of the fracture. Only the inscribed area formed by the arbitrary shaped fracture is used for the analysis. To obtain the desired σ of the fracture, the aperture value h is rescaled. h Fox et al. Geotherm Energy (2015) 3:21 Page 6 of 29 Real reservoirs are not static; there will be dynamic changes to the aperture and per- meability fields due to thermal, hydraulic, mechanical, and chemical (THMC) effects. In earlier works, these effects have been incorporated to model a reservoir’s thermal per - formance for both porous media (McDermott et al. 2006; Taron and Elsworth 2009) and discrete fractures (Ghassemi and Zhou 2011; Rawal and Ghassemi 2014). To achieve via- ble treatments, researchers have to specify quantitative deterministic models of THMC effects which themselves require further additional assumptions. The objective of our study was different and solely directed towards understanding how stochastic spatial aperture variations affect reservoir performance—not to incorporate coupled THMC effects. Ghassemi and Zhou (2011) modeled thermo-poroelastic stresses in a discrete fracture; their work showed that the area near the injection well experienced the largest increase in aperture, driven by local cooling that leads to a large temperature contrast compared to the initial reservoir temperature. In their modeling, the aperture near the injec- tion well increased from 0.23 to 0.38 mm over the operational lifetime of the reservoir. Regions in the middle of the reservoir and near the production well experienced less of an increase. Later, Rawal and Ghassemi (2014) considered silica dissolution, precipita- tion, and diffusion into the surrounding matrix along with thermo-poroelastic effects in a discrete fracture. They demonstrated that these effects changed the original constant aperture of the reservoir, with the cooler region near the injection well increasing by 10 % and the region near the production well decreasing by 24 %. The aperture changes in both of these studies are spatially smooth, which has the effect of lessening any strong perturbations to the base flow and the resulting temperature field. On the other hand, the randomly generated aperture fields with a self-affine description lack smoothness and areas of high and low aperture can occur anywhere in the reservoir, like the sam- ple apertures shown in Fig. 1. Spatial aperture fluctuations modeled as self-affine should not be considered a secondary effect with respect to THMC effects. The magnitude of the initial aperture variation is comparable with the magnitude of the changes due to THMC effects. While THMC induced aperture variations lead to systematic variations in pressure gradient and flow rate, the pre-existing aperture variations produce flow channeling. Fluid flow The fluid flow within the fracture aperture can be modeled as a flow of a Newtonian fluid in a thin channel. For fully developed laminar flow, if the aperture varies gradually, the Navier–Stokes equation is reduced to µ∇ v = ∇P − ρ g, (6) where P is the pressure, v is the velocity, ρ is the density of water, and µ is the dynamic viscosity. The conservation of mass can be expressed as � � ∇ · (ρ w v ) = m ˙ δ(r − r ). w l l (7) The summation on the right hand side refers to sources/sinks of mass, w is the local aperture or gap of the channel, and v is the fluid velocity averaged across the gap. The Fox et al. Geotherm Energy (2015) 3:21 Page 7 of 29 (a) (b) (c) (d) (e) Fig. 1 Dimensionless fracture surface temperature for the base case (a) and fractures exhibiting poor (b) and good (d) thermal performance (left hand side) with their respective aperture fields (c and e; right hand side). The aperture field is normalized by the mean aperture value w and the dimensionless aperture standard deviation σ /w and the Hurst coefficient ζ are equal to 0.25 and 0.8, respectively. The wellbore separation value is L/(2R) = 0.4 and the injection and production wells are located at (−0.4,0) and (0.4,0), respectively (marked by a black dot). All figures represent the same production period of 203 days, the time it takes for the base case to reach a dimensionless production temperature of  = 0.8 p,0 sources and sinks of mass in the fracture come from the connection to the wellbores, the injection (source) and production (sink) wells. Poiseuille flow describes pressure driven laminar flow in a channel with a no slip condition applied at the walls of the channel. For Poiseuille flow, the velocity averaged across the gap is �v� =− ∇P − ρ g , w (8) 12µ Fox et al. Geotherm Energy (2015) 3:21 Page 8 of 29 and this velocity is a function of the position r in the two-dimensional plane of the frac- ture. Combining Eqs. (7) and (8) yields ρ w −∇ · [∇P − ρ g] = m ˙ δ(r − r ), w l l (9) 12µ where the gradient operator only applies in the plane of the fracture. Often, Eq. (9) is referred to as the Reynolds equation and is derived using an approxi- mate lubrication analysis. As such, Eq. (9) is valid when Re(w/L ) ≪ 1 and w/L ≪ 1 x,y x,y (Deen 1998). L refers to the characteristic length scale of aperture variations in the x x,y and y direction and may be considered a “wavelength” of aperture variations. Here, the Reynolds number is defined as Re ≡ (ρ U w)/µ , where U is the characteristic velocity. w c c For a mass flow of 40 kg/s, a wellbore spacing L of 400 m, a width of 1 mm, and a char - acteristic velocity of U ∼˙m/(ρ wL), then Re ∼ 100, well below the critical Reynolds c w number for turbulence transition. The lubrication criteria Re(w/L ) ≪ and w/L ≪ 1 x,y x,y will also be satisfied when Re ∼ 100 because L > 1 m because large scale aperture var- x,y iations dominate when the aperture fluctuations are self-affine. Eq.  (9) has been exten - sively used to model fluid flow in fractures (Brown 1987; Zimmerman and Bodvarsson 1996; Méheust and Schmittbuhl 2001; Watanabe et  al. 2008). Because the aperture is raised to the third power, flow modeled using Eq.  (9) is sometimes said to follow the “cubic law”. The pressure field is first solved using Eq.  (9) with the finite element method (FEM). After solving for the pressure, the velocity is obtained from Eq.  (8). The injection and production wells are modeled as circles of radius 0.15  m to avoid singularities in pres- sure. The mass flux was evenly distributed along the circumference of the injection and production well. We specify a no mass flux condition on the boundary of the fracture and in the direction normal to the surface in order to model the fracture having a lim- ited areal extent with no leak off, that is the surrounding rock matrix is assumed to be impermeable. Heat transport The numerical model we have developed is a hybrid finite element and boundary ele - ment method for discretely fractured geothermal reservoirs. The three-dimensional heat conduction in the rock matrix is captured through a boundary element treatment while the fluid flow and thermal convection are solved via the FEM. As stated earlier, the rock matrix surrounding the discrete fractures is assumed to be impermeable with respect to fluid flow. The numerical model relies on 2D discretization of only the fracture surface and not the entire domain. The 2D discretization offers faster numerical run times over a 3D discretization of the entire reservoir and allows one to invest more computational effort on fracture scale details. The approach to use hybrid FEM and boundary integral representations is not new and our developed model was influenced by such works as Zhou et al. (2009), which studied the effect of thermoelastic stresses in a geothermal res - ervoir. Other hybrid numerical and analytical methods have been derived before, Cheng et al. (2001) used their developed methodology to model multi-dimensional heat effects in geothermal reservoirs and McDermott et al. (2009) used their hybrid model to solve Fox et al. Geotherm Energy (2015) 3:21 Page 9 of 29 for the flow, mass, and heat transport in a fracture, taking advantage of not having to dis - cretize a three dimensional domain. The fluid energy balance inside the fracture is ∂T ρ c w + ρ c w �v� · ∇T − k w∇ T = q . w w w w w s (10) ∂t The term q refers to the surface heat flux and c and k are the heat capacity and ther- s w w mal conductivity of water, respectively. By neglecting the transient heat storage in the liq- uid contained by the fracture and thermal conduction in the liquid (Pe ≡ U L/α ≫ 1 ), c w Eq. (10) becomes ρ c w�v� · ∇T = q . w w s (11) Unfortunately, the preceding equation cannot be solved independently because the surface heat flux is unknown and coupled to transient conduction of the rock mass sur - rounding the fracture. The surface heat flux is also related to the temperature of the sys - tem and can be represented by the 3D Green’s function and the intensity of the source. Here, the source is the fracture surface heat flux. For a source located at z and confined by the extent of the fracture, the temperature of the system is expressed as: t ′ ′ ′ q (x , y , t ) ′ ′ ′ ′ ′ ′ T (x, y, z, t) = T + G(x − x , y − y , z, t − t ) dx dy dt , r (12) ρ c r r 0 A ′ ′ q¯ (x , y , s) ′ ′ ′ ′ ¯ ¯ ¯ T (x, y, z, s) = T + G(x − x , y − y , z, s) dx dy , (13) ρ c r r where G is the 3D Green’s function for transient conduction and the overbar notation refers to the Laplace transform with s representing the Laplace space variable. The sub - script r refers to the thermophysical properties of the rock while T is the original res- ervoir rock temperature. The Laplace transform of the fracture heat balance equation is simply � � ρ c w v · ∇T =¯q . (14) w w s By working in the Laplace space, time discretization is explicitly avoided along with singularities at t = 0 and accuracy issues associated with early time discretization. Equa- tions (13) and (14) are solved simultaneously to obtain the surface heat flux and the tem - perature profile of the fracture at different time steps of interest. The inverse Laplace transform is computed using the Gaver–Stehfest algorithm (Stehfest 1970a, b) which determines the value of f(t) by computing several terms in a series involving f (s) at dif- ferent values of s. ′ ′ The matrix G (x − x , y − y , z) is computationally intensive to construct because every ij entry must be evaluated. The computational cost is O(N ) where N is the number of nodes. Computations can be greatly reduced by noticing that many of the entries in G ij ¯ ¯ are virtually zero and can be skipped in populating G . For determining G for a given ij ij s, entries that correspond to evaluating points that are farther than 20 α/s are skipped where α refers to the thermal diffusivity of the rock. The computational intensity can then be reduced to O(N log N ). The value 20 is arbitrary and can be changed to modify Fox et al. Geotherm Energy (2015) 3:21 Page 10 of 29 accuracy. This method is akin to a single term Fast Multipole Method as described by Greengard and Rokhlin (1997). For a 2500 node mesh, the computational time was reduced by a factor of about 20 with only a difference in the fifth significant figure for the outlet temperature. Alternatively, since the conduction is dominated by the direction perpendicular to the fracture, the 1D Green’s function can be used instead of the full 3D form. Doing so results in a less than 1 % error and over a 20-fold decrease in simulation time compared to using the truncation approximation of G . Additionally, the 1D simplification corre - ij sponds to a 400-fold decrease in simulation time compared to the fully populated G ij matrix. Fox et al. (2013) showed that 1D conduction is a reasonable approximation when β ≪ 1 β . Here, is a dimensionless parameter describing the ratio of the extent of ther- mal depletion normal to the plane of the fracture to the extent of depletion in the plane of the fracture. For dipole flow, β can be defined as β = 2k L/(mc ˙ ), where L is the spac- r w ing between the injection and production well. For typical flow conditions of 40 kg/s and 2 −4 a wellbore separation of 400 m, β = 2 × 10 . Nonetheless, for completeness, we used the full 3D heat conduction in this study. The dimensionless temperature in the reservoir is defined as � ≡ (T(x, y) − T ) /(T − T ), where T is the temperature of the fracture surface at point (x, y), T is the r w w temperature of the injected water and T is the initial rock temperature which remains constant in the far-field. The value of  ranges from zero to one. A value of zero refers to complete thermal exhaustion of the reservoir as the temperature in the fracture is now equal to the injected water temperature. If  = 1, there is no drawdown since the temperature in the fracture is still at the initial reservoir rock temperature. The variable refers to the dimensionless production temperature. If β ≪ 1 (1D conduction), then the dimensionless temperature is a function of the dimensionless wellbore spacing and dimensionless time: L t � = f , , (15) 2R τ 4 2 where τ ≡ L k ρ c /(mc ˙ ) . The dimensionless wellbore spacing is constrained to range r r r w from zero to one. The results for one mass flow rate and a set of times for a given dimen - sionless wellbore spacing are the same as for a different mass flow rate and a set of extrac - tion times corresponding to the same dimensionless times. By adding self-affine aperture fluctuations, the statistical values of the dimensionless temperature become functions of two additional parameters, σ /w and ζ. For this study, we have used a radius R of 500 m and an injected mass flow m ˙ of 40 kg/s. The thermophysical property values used were k = 2.4 W/(m °C), ρ = 2300  kg/m , c = 1000 J/(kg °C), and c = 4184 J/(kg °C). The r r r w value of τ is 2.1 × 10  s or 0.16 years for a wellbore spacing of 400 m and for the previ- ously listed values of the thermophysical properties and mass flow rate. Results and discussion Results for one fracture For an idealized single circular fracture with self-affine aperture variations, our mode - ling provides the statistical distribution of the dimensionless production temperatures as a function of wellbore spacing, duration of extraction, the variance of fracture aperture Fox et al. Geotherm Energy (2015) 3:21 Page 11 of 29 variations, and the Hurst coefficient. The extraction time for geothermal production selected for our study corresponded to a particular level of the thermal performance. In our study, this occurred when the dimensionless production temperature reached 0.8 for the case without aperture variations, which we call the base case. The variable  refers p,0 to the base case production temperature. The drawdown criteria will vary depending on the end-use of the thermal energy; direct thermal use operations can afford a greater drawdown compared to electricity production. However, using  = 0.8 is a realistic p,0 value for the drawdown criteria in many practical situations. For σ /w, ζ, and L/(2R), we chose a range of possible values for these parameters to model their impacts on thermal performance. The range of σ /w was chosen to be σ /w from 0.015 to 0.35. Granite core samples have been shown to have a value of of −4 0.42 (Hakami and Larsson 1996) but also low values around 5 × 10 (Sausse 2002). Because of the uncertainty in the standard deviation of the local aperture values and the difficulty in concluding how representative core sample values will be of practical in situ conditions in geothermal reservoirs, we believe the range chosen is adequate for simulating a large class of fractures. Above σ /�w� = 0.35, it becomes difficult to gener - ate fractures that do not have zero or negative aperture values. u Th s, the aperture reali - zations generated represent fractures that have been sufficiently propped open due to the combination of “self-propping” due to shear displacement and inflation as a result of higher fluid pressures within the fracture. As previously reported, the Hurst coeffi - cient has been observed to be always above 0.5 and is often reported to be around 0.8. However, the Hurst coefficient is a model specific parameter and is not often reported. It is plausible to think that there are instances where the Hurst coefficient is less than 0.5. Given its theoretical range of 0 to 1, the range of the Hurst coefficient chosen was from 0.1 to 1. The dimensionless wellbore spacing L/(2R) can vary anywhere from 0 to 1, depending on where the injection and production wells have been placed relative to the fracture size. In this study, we investigated the range from 0.2 to 0.8. We first present the results where the standard deviation of the aperture is varied. For the range of values investigated, five thousand fracture realizations were generated for each value of σ /w with L/(2R) and ζ kept constant at 0.4 and 0.8, respectively. In a second set of simulations, we repeated the process but now varying the Hurst coeffi - cient and keeping L/(2R) = 0.4 and σ /�w� = 0.25. Finally, in a third set, the wellbore spacing was varied and σ /w and ζ are set to 0.25, and 0.8, respectively. The statistical values of the thermal performance of most interest are the mean dimensionless produc- tion temperature  of the ensemble, the standard deviation of production values σ , p � and the fraction of cases with a higher production temperature than the base case F. Five-thousand was chosen to have an adequate sample size but also to be able to gener- ate results in a timely manner. Using ten-thousand fracture realizations did not result in any significant change in the derived results. Varying strength of aperture fluctuations σ /w For a range of from 0.015 to 0.35, five thousand fracture realizations were gener - ated for each value of σ /w. For each five thousand fracture realization data set, the values of  , σ , and F are reported based on the thermal performance for all realiza- p � tions at the time when the dimensionless production temperature reached 0.8 in the Fox et al. Geotherm Energy (2015) 3:21 Page 12 of 29 absence of aperture variations. A depiction of the fracture surface temperature for the base case, thermally under performing, and thermally over performing (compared with the base case) and their corresponding aperture fields is shown in Fig.  1. The dimension - less heat transfer area accessed at a certain time is defined as (1 − �) dA/A = 1 − ��� . u Th s, the mean fracture surface temperature  is a measure of the thermally accessed area of the fracture at a given extraction time. Comparing  for different fracture aper - ture distributions, a fracture with a lower value of  was able to access more area of the reservoir for heat transfer. For dipole flow, the majority of the flow occurs in the region surrounding the injection and production wells with less flow in the periphery of the fracture. The base case mean value for  is 0.816 when  = 0.8. For a fracture with p,0 smaller apertures in the region between the wells, the fluid is pushed to higher aper - ture regions in the periphery and the effective heat transport area increases, �� = 0.787 when the dimensionless production temperature reaches 0.8. As expected, this less channeled, beneficial flow field leads to a higher production temperature and slows the rate of thermal drawdown. On the other hand, when there exists a high aperture region in between the wells, the flow is further channeled from the injection to production well and the effective heat transfer area becomes narrower than the base case; the mean dimensionless fracture surface temperature is 0.854 when the dimensionless production temperature is 0.8. Consequently, the production well cools down rapidly. How the low and high aperture regions are arranged in the fracture is clearly important in determin- ing whether the reservoir will exhibit poor or good thermal performance. Figure 2 shows the distribution of thermal performance when σ /�w� = 0.25. Although a Gaussian curve does not fit the histogram perfectly, it does an adequate job of describ - ing the behavior. A Gaussian distribution of  is expected for a Gaussian aperture fluc - � � tuations with small amplitude variations. Inserting w + h into the governing equations for fluid velocity and fracture surface temperature and considering a regular perturba - 2 3 tion expansion for small ǫ = σ /�w�, one obtains � = � + ǫ� + ǫ � + O(ǫ ) for w, p 0 1 2 where  is the production temperature for the constant aperture case and  is a linear functional of the field h(r)/σ . A linear functional of a Gaussian field is a Gaussian vari - able. Figure 2 also reveals that we sometimes encounter enhanced thermal performance but more often encounter poor performing fractures with a dimensionless production temperature less than the base case. Figures  3 and 4 show how the strength of aper- ture variations affects the mean and standard deviation of the production temperature. As σ /w increases, the mean production temperature decreases quadratically and the standard deviation increases linearly. The trends reveal that as the aperture fluctuations increase, the thermal performance declines on average and is more variable. In effect, Fig.  4 says that predicting the reservoir’s thermal performance becomes more difficult with increasing values of σ / w . The linear behavior of Fig.  4 is expected due to the before mentioned linear behav- ior of  to small variations of h. In particular,  has zero standard deviation, while p 0 has an O(1) standard deviation so that the standard deviation of � ≈ � + ǫ� is 1 p 0 1 proportional to σ . The quadratic relationship between the mean production tempera - ture and the aperture standard deviation is also expected owing to the fact that  must ¯ ¯ ¯ be zero based on linearity. Thus, � − � = (σ /�w�) � , where  is a constant of 0 h 2 2 proportionality. Fox et al. Geotherm Energy (2015) 3:21 Page 13 of 29 Fig. 2 Histogram of the dimensionless production temperature  when the dimensionless wellbore spac- ing L/(2R) is equal 0.4, the dimensionless standard deviation of the aperture σ /w is equal to 0.25, and the Hurst coefficient ζ is equal to 0.8. The results are for the production period when the base case (no aperture variations) reaches  = 0.8. A Gaussian fit results in a mean and standard deviation of 0.783 and 0.072, p,0 respectively Fig. 3 Deviation of the mean dimensionless production temperature from the base case value  −  as p,0 p a function of the dimensionless standard deviation of the aperture σ /w. The Hurst coefficient ζ is equal to 0.8, the wellbore spacing L/(2R) is equal to 0.4, and the base case production temperature  is equal to 0.8 p,0 (reference value of no aperture variations). A quadratic fit is also shown Geothermal operations are sensitive to temperature drawdown and a premature drop in production temperature would reduce the economic attractiveness of the project. Each fracture in the ensemble has a lifetime corresponding to the time when  reaches 0.8. When σ /�w� = 0.1, which resulted in σ = 0.031, the 25 and 75 % quantiles of the distribution of the ratio of reservoir life to the reservoir life of the base case are equal to 0.83 and 1.16, respectively. For σ /�w� = 0.35, corresponding to σ = 0.092, the results h � are more drastic as the 25 and 75  % quantiles become 0.47 and 1.28, respectively. On the other hand, aperture variations are less important for lower values of σ /w; when � � σ / w = 0.05, the 25 and 75 % quantiles of reservoir life are closer together and equal to 0.91 and 1.08, respectively. Fox et al. Geotherm Energy (2015) 3:21 Page 14 of 29 Fig. 4 Standard deviation of the dimensionless production temperature σ for the ensemble of fracture realizations as a function of dimensionless standard deviation of the aperture σ /w when ζ = 0.8 and L/(2R) = 0.4. The production end time of interest is when the dimensionless production temperature  of the base case (no aperture variations) reaches 0.8. A linear fit is also shown For a normally distributed production temperature, the probability density function is equal to 1 (� − � ) p p f (� ) = exp . p (16) 2σ 2πσ From Fig. 2, it is reasonable to model the production temperature using a normal (Gauss- ian) distribution. The fraction of fractures that exhibit enhanced thermal performance is equal to one minus the integral of Eq. (16) with respect to  form 0 to  . The result is p p,0 1 � − � p,0 p F = erfc √ . (17) 2σ As discussed earlier, the mean and standard deviation of the production temperature distribution have a quadratic and a linear behavior, respectively: � − � = C , (18) p,0 p 1 � � and σ = C . � 2 (19) �w� Using Eqs. (18) and (19), the fraction of fractures that exhibit enhanced thermal per- formance compared to the base case can be expressed in terms of the standard deviation of the aperture fluctuations, 1 Kσ F = erfc √ , (20) 2�w� Fox et al. Geotherm Energy (2015) 3:21 Page 15 of 29 where K is the ratio of C /C . Since the error function is linear for small values of 1 2 K σ /( 2�w�), Eq. (20) can be approximated as 1 K σ F ≈ − √ . (21) 2 �w� 2π Figure  5 displays the linear decrease of F as a function of σ /w. The linear behavior of F results because the argument of the error function is small for the range studied (K = 1.03). The theoretical prediction for F matches well with the simulation results and is consistent with assuming that the distribution of  is Gaussian. The theoretical and simulated curves only deviate about 0.015 at the end of the range investigated. The frac - tion of fractures exhibiting adverse flow channeling was not observed to be a function of dimensionless time as F only varied temporally by less than 1 %. Thus, choosing a differ - ent dimensionless temperature stopping criteria will not change the reported values of F. The results shown indicate that the magnitude of aperture variations has a large effect on the thermal performance of the reservoir. For a specified standard deviation, the pro - duction temperature can be modeled as a Gaussian distribution. The fraction of frac - ture whose thermal performance benefits from a self-affine aperture field monotonically decreases with respect to the strength of the aperture variations and is never higher than 0.5. It ranged from 0.49 to 0.34 for the range of σ /w investigated. Thus, there are more ways to arrange isotropic self-affine aperture fluctuations that lead to a flow field which produces adverse rather than beneficial thermal performance. Varying the Hurst coefficient To investigate the effect of the Hurst coefficient ζ on thermal performance, simulations were performed for ten sets of five thousand fractures. Each set had a different value of ζ in the range from 0.1 to 1 and σ /w was kept constant at 0.25. Figure  6 displays 1D fracture aperture profiles for various values of ζ. Increasing ζ yields aperture fields with strong large scale correlations. In other words, the profile for ζ = 1 is smoother than the Fig. 5 The fraction of fractures exhibiting thermal performance enhancement F as a function of dimen- sionless standard deviation of the aperture σ /w. The Hurst coefficient ζ is equal to 0.8 and the wellbore spacing L/(2R) is equal to 0.4. A best fit line (dashed red) was added as well as the theoretical prediction (solid black) using a Gaussian distribution and the fitted parameters for the relationship of σ and � p Fox et al. Geotherm Energy (2015) 3:21 Page 16 of 29 Fig. 6 The fracture aperture h profile for various values of the Hurst coefficient ζ with the dimensionless standard deviation of aperture fluctuations σ /w kept at 0.25. For values of ζ of 0.1, 0.5, 0.8, and 1, the pro- files were shifted upward by w, 2w, and 3w to prevent them from overlapping profile for ζ = 0.1, which is characterized by many “peaks” and “valleys”. Table  1 sum- marizes the statistical results of the ensembles for each value of ζ, with a 95 % confidence interval derived by resampling based on the bootstrapping method (Efron and Tibshi- rani 1994). Table 1 indicates that the Hurst coefficient has little effect on the thermal performance given the confidence intervals of the statistical quantities. Zimmerman and Bodvarsson (1996) and Méheust and Schmittbuhl (2001) have shown that changing the Hurst coef- ficient had little affect on fluid flow. Since heat transport is mathematically an integral process that is dependent on large length scale fluid flow variations and even the vari - ance of the aperture is governed by the large scale, it is not surprising that varying the Hurst coefficient has little impact on the distribution of thermal performance. Of the statistical values from the ensemble, only σ changes appreciably, ranging from 0.063 to 0.073 with increasing ζ. If the aperture field is more correlated (larger value of ζ), the Table 1 Thermal performance with varying the Hurst coefficient ¯ F ζ σ 0.1 0.786 ± 0.002 0.063 ± 0.001 0.41 ± 0.01 0.2 0.785 ± 0.002 0.065 ± 0.001 0.40 ± 0.01 0.3 0.784 ± 0.002 0.067 ± 0.001 0.40 ± 0.01 0.4 0.784 ± 0.002 0.069 ± 0.001 0.40 ± 0.01 0.5 0.784 ± 0.002 0.069 ± 0.001 0.39 ± 0.01 0.6 0.784 ± 0.002 0.071 ± 0.001 0.40 ± 0.01 0.7 0.783 ± 0.002 0.072 ± 0.001 0.39 ± 0.01 0.8 0.783 ± 0.002 0.072 ± 0.001 0.39 ± 0.01 0.9 0.784 ± 0.002 0.073 ± 0.001 0.40 ± 0.01 1 0.784 ± 0.002 0.073 ± 0.001 0.40 ± 0.01 The statistical values are reported with a 95 % confidence interval determined using resampling based on the bootstrapping method Fox et al. Geotherm Energy (2015) 3:21 Page 17 of 29 region between the wells is more likely to have either a single peak or a single valley (low or high aperture regions). Having either a single peak or valley results in a more drastic difference in thermal performance than would be induced by the many peaks and valleys occurring for smaller ζ. The standard deviation of the production temperature does not change much for ζ> 0.6 because for these values the diameters of the peaks and valleys are likely to be equal or larger than the wellbore spacing. Varying wellbore separation The effect of self-affine aperture variations on different base flows can be understood by changing the value of the wellbore separation value, L/(2R). As L/(2R) gets smaller the streamlines becomes more varied and diverse. On the other hand, the streamlines become more uniform as the injection and production wells are pulled apart. For values of L/(2R) of 0.2, 0.4, 0.5, 0.6, and 0.8, the thermal performance of the five thousand frac - ture simulations for each value of L/(2R) were all performed for an end time when the respective L/(2R) base case production temperature reaches 0.8. To reach the  = 0.8 , the dimensionless time t/τ for L/(2R) of 0.2, 0.4, 0.5, 0.6, and 0.8 is equal to 4.7, 3.5, 2.8, 2.2, and 1.2 respectively. Although increasing the wellbore spacing increases the time to reach the drawdown criteria, the dimensionless time t/τ to reach the drawdown criteria drops for increasing values of wellbore spacing because τ is proportional to L . The val - ues of ζ and σ /w were kept constant at 0.8 and 0.25, respectively. Boxplots for the distribution of thermal performance for the different values of L/(2R) are plotted in Fig.  7. The boxplots show that the distribution of dimensionless produc - tion temperatures becomes more varied with increasing values of the wellbore spacing but peaks at L/(2R) = 0.4; the standard deviation of the dimensionless production tem- perature is greatest for this value. The mean value of the distribution and F decreases as L/(2R) decreases. Table 2 summarizes the statistical values of the data sets, with a 95 % confidence interval derived by resampling based on the bootstrapping method (Efron and Tibshirani 1994). The mean dimensionless production temperature is closest to the base case value when the wellbore spacing is small and there is less opportunity for the conductance to vary in the inter-wellbore region. The thermal performance distribution is most varied when L/(2R) = 0.4 because it corresponds to the typical length scale of one “peak” or “valley”, as can be see in Fig. 6. When the dimensionless wellbore spacing is smaller than 0.4, then the region in between the wells will have less aperture variations which results in a lower σ . The trend of lower values of F for larger values of wellbore spacing can be rationalized by considering the heat transfer area. The mean fracture surface tempera - ture  is a measure of the accessed area for heat transfer at a point in time; smaller val- ues of  mean more area has been accessed. The base case values for  are 0.95, 0.82, 0.74, 0.67 and 0.57 for L/(2R) = 0.2, 0.4, 0.5, 0.6, and 0.8, respectively. For larger wellbore spacing values, the base case is already harvesting energy from a large portion of the fracture and any deflection of the base flow caused by aperture variations will tend to decrease the heat transfer area as seen in the results for F, the fraction of fractures with better performance than the base case. In contrast, smaller wellbore spacing values, the base case heat transfer area is small and perturbations of the flow field are more likely to increase the heat transfer area. Fox et al. Geotherm Energy (2015) 3:21 Page 18 of 29 Fig. 7 Boxplot of dimensionless production temperatures for various values of wellbore spacing L/(2R). The dimensionless standard deviation σ /w is equal to 0.25 and the Hurst coefficient ζ is equal to 0.8. The lower and upper edge of the box represents the first and third quartiles while the solid red line represents the median value. The lower and upper whiskers of the boxplot are defined as 1.5 times the first and third quartile value, respectively, and the red plus signs represent the outliers. The dashed red line represents the base case (no aperture variations) value of 0.8 Table 2 Thermal performance with varying wellbore spacing L/(2R) L/(2R) σ F 0.2 0.793 ± 0.002 0.057 ± 0.001 0.42 ± 0.01 0.4 0.783 ± 0.002 0.072 ± 0.001 0.39 ± 0.01 0.5 0.779 ± 0.002 0.073 ± 0.001 0.38 ± 0.01 0.6 0.774 ± 0.002 0.068 ± 0.001 0.36 ± 0.01 0.8 0.763 ± 0.002 0.053 ± 0.001 0.28 ± 0.01 The statistical values are reported with a 95 % confidence interval determined using resampling based on the bootstrapping method Results for multiple fractures The previous section described the effect of a self-affine aperture field on the thermal performance when there is only one fracture. A relevant situation to further model would be a reservoir consisting of two fractures. We consider a reservoir consisting of two fractures with equal radii that do not intersect with one another, but rather are hydraulically connected via the injection and production well. The pressure drop along the wellbore is assumed negligible because w/r ≪ 1, where the wellbore radius r is around 10 cm and the aperture is around 1 mm. The mass flow was doubled to 80 kg/s, but the value of m ˙ /A remains the same as that in the one fracture system because the total surface area of the system is twice as large. With two fractures the flow field exhibits both inter- and intrafracture flow variations. As has been seen in the one fracture reservoir, local aperture variations disturb the base- flow creating intrafracture flow variations. However, aperture variations can also modify the ratio of the mean flow rate between fractures through altering the flow impedance or resistance to mean flow of each fracture and affect the reservoir’s overall pressure drop. The flow impedance of the fractures in the reservoir is what drives the global distribu - tion of flow in the reservoir. Fox et al. Geotherm Energy (2015) 3:21 Page 19 of 29 Since the pressure drop P is linear with the flow rate, P ∝˙m, the flow split fractions for fracture A and B are φ = , A (22) I + I A B φ = 1 − φ , B A (23) where I ≡ �P/m ˙ and represents the flow impedance of fracture i. The system can be i i viewed as two “resistors” in parallel where the pressure drop is the same between the two fractures. The above equation demonstrates how fracture aperture variations mod - ify the interfracture flow. Since it does not matter whether a fracture is labeled fracture A or B, the labels used in Eqs. (22) and (23) are arbitrary and we may define the variable φ as max(φ , φ ), which ranges from 0.5 to 1. A B The hydraulic aperture is another way to view the varying pressure drop for each frac - ture. In terms of the pressure drop of a fracture P, the hydraulic aperture is defined as 1/3 �P w ≡ w , (24) H 0 �P where P is the reference pressure drop in a smooth fracture for an aperture value of w receiving the same mass flux. The cube root is a result of the cubic law described in Eq. (9). The hydraulic aperture is the aperture value needed in a smooth fracture to have the same pressure drop in a fracture with aperture variations. The hydraulic aperture is a useful metric in quantifying how spatial aperture variations control the pressure drop and affect interfracture flow variations. The effects of fracture aperture variations in two and three fracture systems was stud - ied by determining the thermal performance of an ensemble of fracture realizations. The wellbore spacing and the Hurst coefficient were kept constant and the standard devia - tion of the aperture was varied over the same range considered before. Finally, we exam- ined the efficacy of flow control methods in minimizing any negative effects aperture variations have on the thermal performance of multifracture systems. Varying strength of aperture fluctuations Two thousand five hundred fracture pairs with L/(2R) = 0.4 and ζ = 0.8 were generated and their thermal performance was simulated. Reassigning the pairing had an insignifi - cant effect on the results presented. Figures  8 and 9 show the histogram of the produc- tion temperature and the flow split fraction when σ /�w� = 0.25. It is apparent that the flow split can be modeled as a Gaussian distribution and compared to Fig.  2, the thermal performance is closer to Gaussian than that of the one fracture system. Figure  10 shows that larger φ values, corresponding to more unequal flow splits, are likely to result in lower production temperatures. Also, plotted in the figure is the pro - duction temperature as a function of the flow split if the fractures had no surface vari - ations. In the absence of aperture variations,  is always monotonically decreasing with respect to φ. There are simulation outcomes that lie above the smooth case curve, but systems that exhibit a boost in thermal performance for larger flow split fractions Fox et al. Geotherm Energy (2015) 3:21 Page 20 of 29 Fig. 8 Histogram of the thermal performance of the two fracture system when L/(2R) = 0.4, σ /�w� = 0.25, and ζ = 0.8. A Gaussian fit results in a mean and standard deviation of 0.751 and 0.0625, respectively � � Fig. 9 Histogram of the fractional flow split of the two fracture system when L/(2R) = 0.4, σ / w = 0.25, and ζ = 0.8. By definition the mean is 0.5 and standard deviation is 0.156, respectively become increasingly rare. In general, having a flow split deviate from equality is likely to result in undesirable thermal performance. As plotted in Fig. 11, the uncertainty in ther- mal performance has a linear relationship with the uncertainty in φ. As with the single fracture reservoir, the fraction of reservoirs that benefit from aperture variations decreases with increasing mean square of aperture variations. The fraction of reservoirs with a production temperature higher than the base case, F, as a function of the dimensionless standard deviation, is plotted in Fig.  12. For comparison, the single fracture results are plotted in the same figure. The theoretical fit was obtained by using Eq. (20) along with the fits of and . Compared to the single fracture case, � p the values of both σ and  are smaller. The curve is no longer linear since K is now � p equal to 3.00 and as a result, the linear approximation of the error function breaks down. Fox et al. Geotherm Energy (2015) 3:21 Page 21 of 29 Fig. 10 Dimensionless production temperature as a function of fractional flow split when L/(2R) = 0.4, � � σ / w = 0.25, and ζ = 0.8. The solid black line represents the production temperature for smooth fractures when the flow split varies. One and two standard deviations from the mean are plotted to help visualize how the data is distributed Fig. 11 The standard deviation of the dimensionless production temperature σ as a function of the stand- ard deviation of the fraction flow split σ when L/(2R) = 0.4 and ζ = 0.8. A linear fit is also shown A portion of the decreased thermal performance in the two fracture reservoir results solely from the existence of an unequal flow split between the fractures. The reduction in the production temperature for two smooth fractures with increasing flow split ratio is indicated by the solid line in Fig. 10 and we will refer to this function as g(φ). A metric of the additional loss of thermal performance in the two fracture reservoir compared with the case of two smooth fractures with unequal apertures is F , which we define as the fraction of two fracture reservoirs with � > g(φ), i.e., the fraction that beat the smooth fracture reservoir with the same flow split. The results for F are higher than F because the smooth case thermal performance drops with respect to increasing φ. How- ever, F is closer to the results of F for the two fracture results than to the single fracture results. These results demonstrate that the primary determinant of the diminished per - formance of the two fracture system is the aperture variation and flow channeling within Fox et al. Geotherm Energy (2015) 3:21 Page 22 of 29 Fig. 12 The fraction of two fracture system exhibiting thermal performance enhancement as a function of σ /w. The Hurst coefficient ζ is equal to 0.8 and the wellbore spacing L/(2R) is equal to 0.4. The curves give theoretical predictions based on a Gaussian distribution of  . The squares are the fraction of reservoirs that outperform a base case with an equal flow split, F. The diamonds are the fraction of reservoirs that outper- form a base case with a constant aperture within each fracture but the same flow split as observed in the simulations, F . The one fracture results and the theoretical fit was also plotted each reservoir rather than being simply a result of the difference of hydraulic aperture of the two fractures. The reduced thermal performance compared to a one fracture system can be explained by the relationship between pressure drop and the thermal performance in an individ- ual fracture. In general, a fracture with a small aperture region in between the injec- tion and production wells will exhibit enhanced thermal performance because the flow will expand and sweep through the reservoir more effectively. The reverse is true; a high aperture region in between the wells will promote flow short circuiting and the flow will primarily harvest thermal energy in a narrow area. This behavior was illustrated in Fig.  1, where the temperature of the fracture surface was plotted along side their respective aperture fields. Furthermore, a low aperture wellbore region requires a higher pressure drop to push the fluid through the constricted region surrounding both the injection and production wells. The covariance of w and  for the single fracture is −0.65 when H p L/(2R) = 0.4, ζ = 0.8, and σ /�w� = 0.25, which supports our qualitative argument. For single fracture reservoirs, poor performing fractures exhibit a smaller pressure drop. However, in a two fracture reservoir, the poor thermally performing reservoir will receive a larger portion of the flow and the reservoir with good thermal performance receives a smaller portion of the flow. Thus, it is not simply the existence of an unequal flow split but the fact that the higher flow goes preferentially to the fracture with more channeling that accounts for the poorer performance of two fracture reservoirs com- pared with single fracture reservoirs. However, the negative impact on thermal perfor- mance can be mitigated if the flow split between the reservoirs can be manipulated. Three fracture system We can extend the analysis to a three fracture system, following the same method used for the two fracture reservoir. The system is treated as three resistors in parallel, Fox et al. Geotherm Energy (2015) 3:21 Page 23 of 29 all with the same pressure drop. Plotted in Fig.  13 are the fitted Gaussian probability distribution of the production temperature for one, two, and three fracture reservoir systems for when L/(2R) = 0.4, σ /�w� = 0.25, and ζ = 0.8. With the addition of more fractures, the distribution of thermal performance becomes less varied and the mean shifts to lower values. As expected, the thermal performance of the system gets worse for increasing values of the mean square of aperture fluctuations. For the range of σ /�w� = 0.015 − 0.35 , F monotonically decreases from 0.49 to 0.074. The three frac - ture system has the worst thermal performance among the cases investigated. Using the same argument as the two fracture system, it is more likely that a low impedance frac- ture has individually poor thermal performance and diverts most of the fluid into that fracture. With three fractures, the probability of having the aforementioned outcome increases and thermal performance suffers to a greater extent than the one or two frac - ture reservoirs. Flow control To mitigate the adverse impact of fracture surface variations on the interfracture flow variations, flow control methods can modify the flow split and reduce the negative effects associated with flow split disparity. In a real reservoir, flow control is achieved by using interflow control devices (ICD). An ICD works by increasing the resistance of flow from the wellbore to the formation or fracture. The increased resistance is achieved using nozzles, orifices, and helical channels on the base pipe. For more information on ICD, please refer to Al-Khelaiwi et al. (2010) and Birchenko et al. (2010). Figure  14 shows the boxplot of the dimensionless production temperature distribu- tion of the two fracture system when L/(2R) = 0.4, ζ = 0.8, and σ /�w� = 0.25 for three cases: no flow control (untreated), 50/50 flow split, and the optimal flow split. A 50/50 flow split was chosen as a good guess on how to best distribute the flow. The optimal Fig. 13 Fitted Gaussian probability density function of dimensionless production temperature for one, two, and three fracture system when L/(2R) = 0.4, σ /�w� = 0.25, ζ = 0.8, and with the base case value of = 0.8. The Gaussian distributions were chosen to match the mean and standard deviation of the simu- p,0 lated temperature results Fox et al. Geotherm Energy (2015) 3:21 Page 24 of 29 Fig. 14 Boxplot of dimensionless production temperatures for three cases: no flow control (untreated), 50/50 flow split, and the optimal flow split. The dimensionless standard deviation of the aperture σ /w is kept at 0.25, the Hurst coefficient ζ is equal to 0.8, and the wellbore spacing L/(2R) is equal to 0.4. The lower and upper edge of the box represents the first and third quartiles while the solid red line represents the median value. The lower and upper whiskers of the boxplot are defined as 1.5 times the first and third quartile value, respec- tively, and the red plus sign represents the outliers flow split represents the value of φ that optimizes the production temperature. In prac- tice, the optimal flow split would be very difficult to predict, as one would first need to know the structure of the reservoir from solving an inverse problem based on tracer tests. Nonetheless, it is calculated to be used as an upper bound on thermal perfor- mance. The boxplot of the 50/50 and the optimal flow scenario show that there is only a marginal increase in thermal performance when the optimal value is used. The 25, 50 % (median), and 75 % quantiles of the distribution for the change in production tempera- ture for the equal split are 0.0023, 0.0149, and 0.0458, respectively. For the optimal split, the 25, 50  % (median), and 75  % quantiles of the distribution are 0.0043, 0.0194, and 0.0554, respectively. Another way to view the flow control results is to plot the change in dimensionless production temperature  versus the flow split fraction. For the 50/50 flow split sce - nario, the results for  versus φ are plotted in Fig.  15. The equal flow split strategy does result in worse thermal performance 13 % of the time compared to no flow control. However, 75  % of these cases only saw a minimal drop in the production temperature, less than 0.004. Additionally, using the optimal flow split in these cases led to an increase in  of no more than 0.0019 with improvement occurring in only 75 % of the cases. The figure reveals that as the flow split fraction moves away from a 50/50 split, the improve - ments of applying flow control increases; the covariance of φ and  is 0.88. When the fracture surface is less varied, the distribution of thermal performance is also less varied. Consequently, flow control will become less important. When σ /�w� = 0.05, 75  % of the cases only had a  of less than 0.002 when the optimal value of φ was used. Compared to the F value for the one fracture reservoir, applying flow control for the range of σ /w investigated was capable of increasing F to a range of 0.50–0.31, which represents 100–89 % of the value of the one fracture reservoir. Plotted in Fig. 16 are the F values for the 50/50 flow split as well as the values for the one and two fracture cases (no flow control). It is only when σ /w > 0.25 that F of the equal flow split is appreciably lower than the one fracture reservoir. Fox et al. Geotherm Energy (2015) 3:21 Page 25 of 29 Fig. 15 The change of production temperature  , using an equal flow split strategy as a function of the original flow split fraction when L/(2R) = 0.4, σ /�w� = 0.25, and ζ = 0.8. The original flow split is normally distributed and one and two standard deviations from the mean are plotted to help visualize how  is distributed The fact that flow control was able to bring F near the single fracture results indicates that interfracture flow variations inhibit thermal performance more than intrafracture flow variations. The results in Fig.  12 showed that these detrimental interfracture flow variations consist of flow preferentially going to poor thermal performing fractures rather than simply involving unequal flow splits. Nonetheless a comparison of our two flow control strategies indicates that an effort to divert flow preferentially away from fractures exhibiting channeling yields only a marginal improvement over a flow control strategy that yields an equal flow split. The flow control results reveal that for equally sized fractures, a 50/50 flow split is good practice. Knowing the optimal split does not result in a large improvement in ther- mal performance. Flow control would be limited if the reservoir consists of fracture- to-fracture intersections as ICD can only control flow in the fractures that are directly Fig. 16 The fraction of two fracture system exhibiting thermal performance enhancement as a function of σ /w for the case of equal flow split. The Hurst coefficient ζ is equal to 0.8 and the wellbore spacing L/(2R) is equal to 0.4. The results for the one fracture, and untreated two fracture cases are also plotted Fox et al. Geotherm Energy (2015) 3:21 Page 26 of 29 connected to the wellbore. Fracture-to-fracture intersections, therefore, open up new possible flow paths rather than just partitioning the flow between two independent flow paths where the path with the least impedance takes the majority of the flow. Thus, the presence of fracture-to-fracture intersections may not result in a diminished thermal performance, as was shown in the case of multiple fractures only hydraulically con- nected via the wellbores. Conclusions Isotropic self-affine fracture aperture variations limit a reservoir’s thermal performance. There are more ways to arrange aperture variations that causes flow paths that lead to premature drawdown than flow pathways that enhance heat transfer. We have shown σ /w that the dimensionless standard deviation of the aperture is more important in controlling thermal performance than the Hurst coefficient ζ. Specifically, both the uncertainty in thermal performance, σ , and the number of fracture realizations that benefit from aperture variations have linear relationships to the root mean square of the aperture variations. However, the effects of aperture variations become nearly insignifi - cant when σ /w < 0.1. For the investigated range of σ /w from 0.015 to 0.35, 34 to h h 49  % of the fracture realizations benefited from aperture variations. Reservoirs with a larger wellbore separation have more uniform base flow, flow in the absence of aper - ture variations, that is more negatively affected by aperture variations. A reservoir with a large wellbore separation already accesses a larger area of the fracture for heat exchange and any flow channeling induced by aperture variations will more likely decrease the effective heat exchange area. Our finding that aperture variations will lead to adverse thermal performance, more often than not, is in agreement with the previous work of Neuville et al. (2010a, b), despite the fact that their study did not couple heat transport in the fracture with that in the surrounding rock and used a uniform base flow in a rectan - gular fracture. Additionally, our use of statistical values for an ensemble of realizations, including the mean and standard deviation of production temperatures, offers a better and more practical understanding of the effects fracture aperture variations have on dis - cretely fractured geothermal reservoirs. We also demonstrated that fracture aperture variations alter the interfracture flow split in a system with multiple parallel fractures. The number of reservoirs exhibiting diminished thermal performance increased compared a reservoir with one fracture and the same mass flow per fracture area. Simply having a higher flow rate in one fracture will more than likely lead to a worse thermal performance because the fracture receiving most of the flow will experience thermal breakthrough more quickly. Further deteriora - tion in thermal performance arises because fractures with a larger flow impedance usu - ally have better thermal performance but receive a smaller portion of the flow. However, applying flow control can alleviate the problems with undesirable flow splits, suggesting a need for inflow control devices specifically suited for geothermal applications rather than “off the shelf ” devices made for oil and natural gas operations. Our analysis suggests that a major cause of poor thermal performing reservoirs is due to interfracture flow variations. With more fractures, the reservoir is more susceptible to adverse thermal performance. In general, more fractures in the reservoir lead to a drop in the variability of thermal performance but also leads to a higher likelihood of having Fox et al. Geotherm Energy (2015) 3:21 Page 27 of 29 a worse production temperature compared to the base case scenario of no aperture variations. The multiple fracture systems simulated in this study were only hydraulically connected via the injection and production well. Real EGS reservoirs, however, like the one at Soultz-sous-Forêts will have fractures that intersect with each other (Sausse et al. 2010). With fracture-fracture intersections, flow control devices become less effective because they can only control well-to-fracture flow. The present study dealt only with a static reservoir where the spatially varying aper - ture field stayed constant with time. Real reservoirs, however, will experience ther - mal, hydraulic, mechanical, and chemical (THMC) changes that alter the aperture field throughout the energy extraction process. In addition to THMC effects, fracture shearing/slipping during heat extraction will also change the aperture field. In general, aperture change due to THMC effects will increase permeability/aperture in the region enclosed by the injection and production well (Taron and Elsworth 2009; Ghassemi and Zhou 2011). Fracture slippage will lead to an increase in the standard deviation of the aperture fluctuations and result in preferential flow channeling in the direction perpen - dicular to the shear displacement (Yeo et  al. 1998; Plouraboué et  al. 2000; Watanabe et al. 2008). Although our work did not incorporate fracture slippage and THMC effects it, nonetheless, provides a framework through which to understand how spatial aperture fluctuations affect thermal performance. It shows that fracture surface variations are a reservoir characteristic that researchers should consider in modeling discretely fractured geothermal reservoirs and performing techo-economic analysis (Beckers et  al. 2014; Held et al. 2014). Our systematic approach for generating spatially varying apertures can be used in conjunction with modeling the effects that THMC and fracture slippage have on reservoir thermal performance, by providing a base aperture field that already has a spatially varying aperture field. The developed numerical framework and the lessons learned from this study are cur - rently being applied to a field laboratory at the Altona Flat Rocks in Plattsburgh, New York. The Altona reservoir represents a meso-scale field site, larger than a bench scale, artificial reservoir but smaller than a full sized geothermal reservoir. The Altona reser - voir consists of Cambrian Potsdam Sandstone formation, where a five spot well system has been drilled in a 10 m by 10 m area, intersecting an isolated horizontal fracture 7.6 m deep (Becker and Tsoflias 2010). Ground penetrating radar has been used to determine the spatial aperture of the fracture (Tsoflias and Becker 2008), creating a reservoir with a characterized subsurface. The site has been used to conduct tracer tests, demonstrating the presence of flow channelizations caused by the spatially varying aperture field (Cast - agna et al. 2011; Becker et al. 2015; Hawkins et al. 2015). The site resembles the reservoir model used in this study, with an induced dipole flow, Gaussian distribution of aper - ture, and a self-affine spatial correlation. Future experimentation at the site will include deployment of thermally degrading and sorbing tracers and reservoir heating through injection of hot water as an analog of the geothermal energy extraction process. Authors’ contributions This study was performed by DBF with advice on the planning, design, and geothermal application of the model pro- vided by DLK and JWT. DBF wrote the manuscript and it was revised and the final version approved by all three authors. All authors read and approved the final manuscript. Author details 1 2 Chemical and Biomolecular Engineering, Cornell University, Ithaca 14850, USA. Cornell Energy Institute, Cornell University, Ithaca 14850, USA. Fox et al. Geotherm Energy (2015) 3:21 Page 28 of 29 Acknowledgements The authors would like to thank the National Science Foundation IGERT program (contract 0966045), the U.S. Depart- ment of Energy (contract DE-EE0000842), and the Cornell Energy Institute for partial financial support for this work. They would also like to thank their colleagues at Cornell University, specifically Koenraad Beckers, Adam Hawkins, Maciej Łukawski, and Calvin Whealton for their input and suggestions about this work. The authors appreciate the feedback and advice received from the two anonymous reviewers. Their input helped strengthen this work. Competing interests The authors declare that they have no competing interests. Received: 5 August 2015 Accepted: 30 September 2015 References Al-Khelaiwi FT, Birchenko VM, Konopczynski MR, Davies D. Advanced wells: a comprehensive approach to the selection between passive and active inflow-control completions. SPE Prod Oper. 2010;25:305–26. Auradou H, Drazer G, Boschan A, Hulin J-P, Koplik J. Flow channeling in a single fracture induced by shear displacement. Geothermics. 2006;35(5–6):576–88. Auradou H, Boschan A, Chertcoff R, DAngelo, M-V, Hulin J-P, Ippolito I. Miscible transfer of solute in different model frac- tures: From random to multiscale wall roughness. Comptes Rendus Geoscience. 2010;342(7–8):644–52. Becker MW, Tsoflias GP. Comparing flux-averaged and resident concentration in a fractured bedrock using ground pen- etrating radar. Water Resour Res. 2010;46(9):1–12 Becker MW, Tsoflias GP, Baker M, Hawkins A. Confirmation of hydraulic, tracer, and heat transfer characterization of a frac- tured bedrock using ground penetrating radar. In: Proceedings of the Fortieth Workshop on Geothermal Reservoir Engineering, Stanford University, January 26–28, 2015, Stanford, CA. 2015. Beckers KF, Lukawski MZ, Anderson BJ, Moore MC, Tester JW. Levelized costs of electricity and direct-use heat from Enhanced Geothermal Systems. J Renew Sustain Energy. 2014;6(1):1–15 Birchenko VM, Muradov KM, Davies DR. Reduction of the horizontal well’s heeltoe effect with inflow control devices. J Pet Sci Eng. 2010;75:244–50. Boffa JM, Allain C, Chertcoff R, Hulin JP, Plouraboué F, Roux S. Roughness of sandstone fracture surfaces: Profilometry and shadow length investigations. Eur Phys J B. 1999;7(2):179–82. Bouchaud E. Scaling properties of cracks. J Phys Condens Matter. 1997;9(21):4319–44. Brown SR. Fluid flow through rock joints: the effect of surface roughness. J Geophys Res Solid Earth. 1987;92(B2):1337–47. Brown SR, Scholz CH. Broad bandwidth study of the topography of natural rock surfaces. J Geoph Res Solid Earth. 1985;90(B14):12575–82. Brown SR, Kranz RL, Bonner BP. Correlation between the surfaces of natural rock joints. Geophy Res Lett. 1986;13(13):1430–3. Bruderer-Weng C, Cowie P, Bernabé Y, Main I. Relating flow channelling to tracer dispersion in heterogeneous networks. Adv Water Resour. 2004;27(8):843–55. Bruel D, Cacas MC, Ledoux E, de Marsily G. Modelling storage behaviour in a fractured rock mass. J Hydrol. 1994;162(3–4):267–78. Cacas MC, Ledoux E, de Marsily G, Tillie B, Barbreau A, Durand E, Feuga B, Peaudecerf P. Modeling fracture flow with a sto - chastic discrete fracture network: Calibration and validation 1. The flow model. Water Resour Res. 1990;26(3):479–89. Castagna M, Becker MW, Bellin A. Joint estimation of transmissivity and storativity in a bedrock fracture. Water Resour Res. 2011;47(9):1–19. Cheng AHD, Ghassemi A, Detournay E. Integral equation solution of heat extraction from a fracture in hot dry rock. Int J Numer Anal Methods Geomech. 2001;25(13):1327–38. Deen WM. Analysis of Transport Phenomena. New York: Oxford University Press; 1998. Drazer G, Koplik J. Permeability of self-affine rough fractures. Phys Rev E. 2000;62(6):8076–85. Drazer G, Koplik J. Transport in rough self-affine fractures. Phys Rev E. 2002;66:026303. Efron B, Tibshirani RJ. An Introduction to the Bootstrap. New York: Chapman & Hall; 1994. Falconer K. Fractal geometry: mathematical foundations and application. 2nd ed. Chichester: Wiley; 2003. Feder J. Fractals physics of solids and liquids. New York: Plenum Press; 1988. Fox DB, Sutter D, Beckers KF, Lukawski MZ, Koch DL, Anderson BJ, Tester JW. Sustainable heat farming: Modeling extrac- tion and recovery in discretely fractured geothermal reservoirs. Geothermics. 2013;46:42–54. Ghassemi A, Zhou X. A three-dimensional thermo-poroelastic model for fracture response to injection/extraction in enhanced geothermal systems. Geothermics. 2011;40(1):39–49. Glover PWJ, Matsuki K, Hikima R, Hayashi K. Fluid flow in synthetic rough fractures and application to the Hachimantai geothermal hot dry rock test site. J Geophys Res Solid Earth. 1998;103(B5):9621–35. Greengard L, Rokhlin V. A fast algorithm for particle simulations. J Comput Phys. 1997;135(2):280–92. Hakami E, Larsson E. Aperture measurements and flow experiments on a single natural fracture. Int J Rock Mech Min Sci Geomech Abstr. 1996;33(4):395–404. Hawkins A, Fox D, Zhao R, Tester J, Cathles L, Koch D, Becker M. Predicting thermal breakthrough from tracer tests: Simulations and observations in a low temperature field laboratory. In: Proceedings of the Fortieth Workshop on Geothermal Reservoir Engineering, Stanford University, January 26–28, 2015, Stanford, CA. 2015. Held S, Genter A, Kohl T, Koelbel T, Sausse J, Schoenball M. Economic evaluation of geothermal reservoir performance through modeling the complexity of the operating EGS in Soultz-sous-Forêts. Geothermics. 2014;51:270–80. Jupe AJ, Bruel D, Hicks T, Hopkirk R, Kappelmeyer O, Kohl T, Kolditz O, Rodrigues N, Smolka K, Willis-Richards J, Wallroth T, Xu S. Modelling of a European prototype HDR reservoir. Geothermics. 1995;24(3):403–19. Fox et al. Geotherm Energy (2015) 3:21 Page 29 of 29 Kolditz O. Modelling flow and heat transfer in fractured rocks: Dimensional effect of matrix heat diffusion. Geothermics. 1995;24(3):421–37. Mandelbrot BB, Ness JWV. Fractional Brownian motions, fractional noises and applications. SIAM Rev. 1968;10(4):422–37. McDermott CI, Randriamanjatosoa ARL, Tenzer H, Kolditz O. Simulation of heat extraction from crystalline rocks: the influ- ence of coupled processes on differential reservoir cooling. Geothermics. 2006;35(3):321–44. McDermott CI, Walsh R, Mettier R, Kosakowski G, Kolditz O. Hybrid analytical and finite element numerical modeling of mass and heat transport in fractured rocks with matrix diffusion. Comput Geosci. 2009;13(3):349–61. Méheust Y, Schmittbuhl J. Geometrical heterogeneities and permeability anisotropy of rough fractures. J Geophys Res Solid Earth. 2001;106(B2):2089–102. Méheust Y, Schmittbuhl J. Scale effects related to flow in rough fractures. Pure Appl Geophys. 2003;160(5–6):1023–1050. Neuville A, Toussaint R, Schmittbuhl J. Fracture roughness and thermal exchange: a case study at Soultz-sous-Forêts. Comptes Rendus Geosci. 2010a;342(7–8):616–25. Neuville A, Toussaint R, Schmittbuhl J. Hydrothermal coupling in a self-affine rough fracture. Phys Rev E. 2010b;82:036317. Nicol DAC, Robinson BA. Modelling the heat extraction from the Rosemanowes HDR reservoir. Geothermics. 1990;19(3):247–57. Nordqvist AW, Tsang YW, Tsang CF, Dverstorp B, Andersson J. A variable aperture fracture network model for flow and transport in fractured rocks. Water Resour Res. 1992;28(6):1703–13. Peitgen HO, Saupe D. The Science of Fractal Images. New York: Springer; 1988. Plouraboué F, Kurowski P, Hulin J-P, Roux S, Schmittbuhl J. Aperture of rough cracks. Phys Rev E. 1995;51(3):1675. Plouraboué F, Hulin J-P, Roux S, Koplik J. Numerical study of geometrical dispersion in self-affine rough fractures. Phys Rev E. 1998;58:3334–46. Plouraboué F, Kurowski P, Boffa J-M, Hulin J-P, Roux S. Experimental study of the transport properties of rough self-affine fractures. J Contam Hydrol. 2000;46(3–4):295–318. Ponson L, Auradou H, Pessel M, Lazarus V, Hulin J-P. Failure mechanisms and surface roughness statistics of fractured fontainebleau sandstone. Phys Rev E. 2007;76(3):036108. Rawal C, Ghassemi A. A reactive thermo-poroelastic analysis of water injection into an enhanced geothermal reservoir. Geothermics. 2014;50:10–23. Sausse J. Hydromechanical properties and alteration of natural fracture surfaces in the Soultz granite (Bas-Rhin, France). Tectonophysics. 2002;348(1–3):169–85. Sausse J, Dezayes C, Dorbath L, Genter A, Place J. 3D model of fracture zones at Soultz-sous-Forêts based on geo- logical data, image logs, induced microseismicity and vertical seismic profiles. Comptes Rendus Geosci. 2010;342(7–8):531–45. Schmittbuhl J, Gentier S, Roux S. Field measurements of the roughness of fault surfaces. Geophys Res Lett. 1993;20(8):639–41. Schmittbuhl J, Schmitt F, Scholz C. Scaling invariance of crack surfaces. J Geophys Res Solid Earth. 1995;100(B4):5953–73. Stehfest H. Algorithm 368: Numerical inversion of Laplace transforms [d5]. Commun ACM. 1970a;13(1):47–9. Stehfest H. Remark on algorithm 368: Numerical inversion of Laplace transforms. Commun ACM. 1970b;13(10):624. Taron J, Elsworth D. Thermal-hydrologic-mechanical-chemical processes in the evolution of engineered geothermal reservoirs. Int J Rock Mech Min Sci. 2009;46(5):855–64. Tsang YW, Tsang CF. Flow channeling in a single fracture as a two-dimensional strongly heterogeneous permeable medium. Water Resour Res. 1989;25(9):2076–80. Tsoflias GP, Becker MW. Ground-penetrating-radar response to fracture-fluid salinity: Why lower frequencies are favorable for resolving salinity changes. Geophysics. 2008;73(5):25–30. Watanabe N, Hirano N, Tsuchiya N. Determination of aperture structure and fluid flow in a rock fracture by high-resolution numeri- cal modeling on the basis of a flow-through experiment under confining pressure. Water Resour Res. 2008;44(6):1–11. Watanabe N, Wang W, McDermott CI, Taniguchi T, Kolditz O. Uncertainty analysis of thermo-hydro-mechanical coupled processes in heterogeneous porous media. Comput Mech. 2010;45(4):263–80. Yeo IW, de Freitas MH, Zimmerman RW. Eec ff t of shear displacement on the aperture and permeability of a rock fracture. Int J Rock Mech Min Sci. 1998;35(8):1051–70. Zhou XX, Ghassemi A, Cheng AH-D. A three-dimensional integral equation model for calculating poro- and thermoe- lastic stresses induced by cold water injection into a geothermal reservoir. Int J Numer Anal Methods Geomech. 2009;33(14):1613–40. Zimmerman RW, Bodvarsson GS. Hydraulic conductivity of rock fractures. Transp Porous Media. 1996;23(1):1–30. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geothermal Energy Springer Journals

The effect of spatial aperture variations on the thermal performance of discretely fractured geothermal reservoirs

Geothermal Energy , Volume 3 (1) – Oct 28, 2015

Loading next page...
 
/lp/springer-journals/the-effect-of-spatial-aperture-variations-on-the-thermal-performance-luBOL3otR7
Publisher
Springer Journals
Copyright
Copyright © 2015 by Fox et al.
Subject
Earth Sciences; Geotechnical Engineering & Applied Earth Sciences; Renewable and Green Energy; Geoecology/Natural Processes
eISSN
2195-9706
DOI
10.1186/s40517-015-0039-z
Publisher site
See Article on Publisher Site

Abstract

cornell.edu Cornell Energy Institute, The effect of spatial aperture variations on the thermal performance of discretely frac- Cornell University, tured geothermal reservoirs was investigated using finite element method solutions to Ithaca 14850, USA the convective heat transfer in the fracture coupled with a boundary integral descrip- Full list of author information is available at the end of the tion of conductive heat transfer in the rock. The dipolar flow was generated by a source article of fluid volume at an injection well and a sink at a production well within a circular fracture. The statistics of the thermal performance of an ensemble of fracture realiza- tions was evaluated. Fractures with self-affine aperture fields where long range correla- tions dominant over short range correlations were generated. The results showed that spatial aperture variations most frequently lead to diminished thermal performance by creating flow channeling that reduces the heat transfer area. Enhanced thermal performance occurred in the less common cases when the aperture was small in the region between the wellbores, causing fluid flow to sweep out greater areas of the fracture and extract heat from a larger area. The standard deviation of the apertures had the largest influence on the thermal performance, while the spatial correlation of the aperture played a secondary role. Larger values of the standard deviation led to more adverse thermal performance. For the range of standard deviations investigated, the fraction of fractures exhibiting enhanced thermal performance compared to the base case of no aperture variations ranged from 34 to 49 %. The degradation of thermal performance due to aperture variations was largest when the well bore spacing was a larger fraction of the fracture diameter. Reservoirs consisting of two non-intersecting fractures connected to the same injection and production wells were also modeled. The uneven split of the flow between the reservoirs in this case caused a further deterioration of thermal performance compared with the single fracture reservoirs. However, flow control was able to overcome nearly all of the additional loss of thermal performance for the multifracture reservoir. Keywords: Enhanced geothermal systems, Aperture variations, Roughness, Discrete fractures, Heat transport in fractures , Coupled convective conductive heat transport, Thermal hydraulic modeling, Geothermal reservoir performance modeling Background Variations in fluid velocity on various flow paths in discretely fractured geothermal res - ervoirs like those found in Enhanced/Engineered Geothermal Systems (EGS) have a sig- nificant effect on the thermal performance of a reservoir. One source of variations in © 2015 Fox et al. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Fox et al. Geotherm Energy (2015) 3:21 Page 2 of 29 fluid flow is from spatial aperture variations in the fractures. EGS reservoirs are devel - oped in hot rock with low permeability and porosity to extract a portion of the stored thermal energy contained in the rock mass. A reservoir is stimulated by pressuriza- tion to open existing sealed fractures (hydro-shearing) or induce new ones by hydraulic fracturing where water may be pumped into and out of the fractured reservoir using a connected set of injection and production wells to “farm” or extract the stored ther- mal energy. While convective heat transport is the dominant mode of transport within the fractures, it is inherently coupled to conductive heat transport from the surrounding −6 2 bulk rock. Because of the low thermal diffusivity of the host rock (of order 10 m /s), only a small region surrounding the fracture network is effective for heat extraction. Given the finite areal extent of the fracture network, extended use will locally cool down the reservoir. The mass loading per unit area usually expressed as the total flow rate per unit fracture surface area is a critical parameter that often scales directly with the rate of thermal drawdown. Drawdown is a major concern for geothermal developers who must forecast reservoir performance and consider its effect in their techno-economic analysis. This study seeks to determine the effect of fracture aperture variations on the thermal performance of discretely fractured geothermal reservoirs Because fractures do not have smooth surfaces, spatial aperture variations can lead to maldistributed flow within the fracture network adding a degree of thermal performance uncertainty. Fracture core samples have revealed fracture apertures to have a self-aff - ine spatial correlation governed by a roughness or Hurst coefficient (Brown et al. 1986; Glover et al. 1998; Schmittbuhl et al. 1993, 1995; Plouraboué et al. 1995; Boffa et al. 1999; Ponson et  al. 2007). A self-affine spatial correlation is one characterized by long range correlations. The aperture undulations in the core samples are found on all length scales but the aperture variance is dominated by the largest scale. Often self-affine fracture surface variations are referred to as fracture roughness despite the fact that roughness implies variations on small rather than large length scales. In this paper, we have avoided using “rough and roughness” as descriptors as much as possible and instead use terms like aperture or surface variations. To determine the effect fracture surface variations have on fluid flow and on a reser - voir’s thermal performance, a dipole flow in circular fractures between a single injector and producer was selected. Dipole flow is not uniform and is characterized by stream - lines where the fluid velocity is faster near the injection and production wells. For sim - plicity, effects like temperature dependent properties and buoyancy were not considered. We are interested in modeling reservoirs with different degrees of fracture surface varia - tions, roughness/Hurst coefficients, wellbore spacings, and multiple fractures. Our focus is to determine how fracture aperture variations affect reservoir parameters like the pro - duction temperature and the likelihood a reservoir will exhibit thermal enhancement or decline through statistical parameters involving an ensemble set of fracture realizations. Through these statistical values, we hope to 1. infer how the extent of aperture variations aid or hinder thermal performance, 2. characterize the degree of reservoir performance uncertainty due to aperture variations, and 3. determine under what scenarios aperture variations need to be modeled. Furthermore, modeling spatial aperture variations pro- vides a method to systematically generate variations in flow fluid rates and their pathway through the fracture system. Fox et al. Geotherm Energy (2015) 3:21 Page 3 of 29 While little work has been done on the coupled convective-conductive heat trans- port of the fracture and surrounding rock in a reservoir with spatial aperture variations, researchers have modeled fluid flow and other transport properties using a self-affine aperture field. They have observed that self-affine aperture fields lead to flow chan - neling where most of the fluid will travel through several preferential channels (Brown 1987; Tsang and Tsang 1989; Glover et al. 1998; Drazer and Koplik 2000; Méheust and Schmittbuhl 2001, 2003; Auradou et  al. 2006). Méheust and Schmittbuhl (2001) deter- mined that three parameters are required to model the flow, the mean aperture, the root mean square of the aperture, and the Hurst coefficient. They further showed how the mean flow rate was either enhanced or inhibited, depending on the orientation of the applied pressure drop. Transport properties such as dispersion have also been studied in self-affine fractures (Plouraboué et al. 1998; Drazer and Koplik 2002; McDermott et al. 2009; Auradou et al. 2010). Numerical studies exhibit non-Gaussian transit time distri- bution with long tails (Drazer and Koplik 2002) and experiments with Newtonian and shear thinning fluids showed that geometric/hydraulic dispersion is dominant at low Péclet numbers with the dispersion coefficient on the order of the Péclet number (Aura - dou et al. 2010). Other researchers have used flow channeling to reduce numerical computation by reducing flow in 2D fractures to an equivalent 1D flow (Cacas et al. 1990; Nordqvist et al. 1992; Bruel et  al. 1994) or to reduce 3D flow into an equivalent 2D system (Bruderer- Weng et al. 2004). The above mentioned studies have focused on the effects fracture sur - face variations have on hydraulic performance and particle/solute transport and do not cover the consequence on convective heat transport. However, these studies mention preferential flow paths which suggest that surface variations may have a negative impact on thermal exchange. Kolditz (1995) suggested that local aperture variations may play a strong role in a res- ervoir’s thermal performance. Jupe et al. (1995) attempted to model spatial aperture var- iations by modeling separate flow paths. They used the multifracture parallel rectangular reservoir model described by Nicol and Robinson (1990) and treated the different flow paths as separate fractures. Their work focused on determining the wellbore spacing needed to obtain an economic EGS reservoir. By incorporating multiple flow paths, the wellbore spacing needed to obtain their desired performance criteria went from 700 m (no channeling) to 1080 m. The uncertainty of reservoir parameters like heat capacity and permeability have been modeled for a heterogeneous porous media geothermal system (Watanabe et al. 2010). The spatial correlation of the reservoir parameters were modeled using a spherical vari - ogram. Using a hundred realizations of the reservoir, Watanabe et al. (2010) determined that permeability was the most sensitivity parameter, with the production temperature having an uncertainty range of 40 K for the 15 years. Since the aperture can be viewed as an analog to permeability, we would expect spatial aperture variations to substantially control a reservoir thermal performance. The thermal effect of surface variations has recently been investigated in a rectangular fracture where the base flow is uniform (Neuville et al. 2010a , b). However, Neuville et al. (2010a, b) used a simplified heat transfer model that only considers the heat transfer in the fracture and not in the surrounding rock; the surrounding rock was assumed to be Fox et al. Geotherm Energy (2015) 3:21 Page 4 of 29 permanently at the original reservoir temperature. As a consequence, their results can be interpreted as the early thermal behavior of the reservoir and not the behavior after extended reservoir use. The assumption of constant rock temperature is valid for the time scale of w /α, where w is the fracture aperture and α is the thermal diffusivity of the −6 2 rock. For an aperture value from 1 to 5 mm and a thermal diffusivity of 1 × 10  m /s, the time scale will be on the order of 1–10  s. Neuville et  al. (2010b) work focused on developing relationships between a reservoir’s hydraulic performance and its thermal performance. They concluded heat exchange in self-affine fractures is less efficient com - pared to a constant aperture case. As a case study, Neuville et  al. (2010a) modeled the Soultz-sous-Forêts EGS reservoir as a single fracture and concluded that merely know- ing the mean aperture is not sufficient to determine the thermal heat exchange and that fracture surface variations always led to inhibited thermal exchange. Although Neuville et al. (2010a, b) have modeled variations in convective heat transport in the presence of aperture undulations, there is still a need to model the coupled convective heat trans- port in the fracture with that of the conductive transport in the surrounding rock. Our approach fully couples the heat transfer in the fracture to the surrounding bulk rock by using a finite element method in the fracture and boundary integral equations to describe the temperature in the bulk rock. Additionally, our study aims to determine how much deviation in production temperature in the absence of aperture undulations is anticipated after prolonged reservoir utilization. The use of statistical values of an ensemble of realizations including the mean and standard deviation of production tem- peratures will offer a better and more practical understanding of the effects fracture sur - face variations have on discretely fractured geothermal reservoirs. Methods The following section outlines our numerical approach in modeling self-affine aperture fields and the thermal performance of the resulting reservoirs. We first introduce the mathematical theory of self-affine aperture fields and describe the algorithm used to generate the fracture surfaces. Then, we detail the method used to solve for the fluid flow in the fracture and finally, we describe our numerical approach in solving for the heat transport in the reservoir and the transient temperature field. Self‑affine aperture field A spatially varying aperture field fluctuates with a value of h(x,  y) from a mean value of w to give w(x, y) = �w� + h(x, y). (1) The spatial mean aperture w is given by w(x, y) dA �w� ≡ , (2) dA where A is the fracture surface area and the variance or mean square of the aperture field is given by h (x, y) dA 2 2 A σ ≡ h = . dA A (3) Fox et al. Geotherm Energy (2015) 3:21 Page 5 of 29 The spatial correlation of an aperture field can be described by the variogram function γ, which is twice the difference of the aperture variance and the two-position autocor - relation function: ′ ′ 2 ′ γ(r ) ≡ h(r) − h(r + r ) = 2 h − h(r)h(r + r ) , (4) where r is the displacement from point r and r is a point on the surface of the fracture. The autocorrelation function is a measure of the correlation of the aperture at point r with that at point r + r . If the variogram function follows ′ ′ 2ζ γ(r ) = C|r | , (5) ′ 2ζ ′ it has a scaling behavior of γ(r ) =  γ(r ), where  is a scaling parameter and the field is said to have a self-affine geometry governed by the Hurst or roughness coefficient, ζ (Feder 1988; Falconer 2003). For this study, we have chosen to model fracture aperture variations as having a self-affine geometry. By measuring a surface profile, many mate - rials have been shown to have a self-affine geometry (Bouchaud 1997). From fracture core samples, the Hurst coefficient has been determined to be from 0.7 to 0.8 (Brown et al. 1986; Glover et al. 1998) and, specifically, equal to 0.8 for granite (Schmittbuhl et al. 1993, 1995; Plouraboué et  al. 1995), the typical rock type of deep geothermal systems. Sandstone, on the other hand, has been shown to have a lower value of ζ at around 0.5 (Boffa et  al. 1999; Ponson et  al. 2007). Although it is difficult to determine whether a fracture will exhibit a self-affine behavior at the length scale of geothermal fractures, it has been observed up to the meter scale (Brown and Scholz 1985). With a self-affine behavior, larger scale variations in h (small frequencies, long wave - lengths) are more important than small scale (large frequencies, short wavelengths) vari- ations. The large scale variations become more strongly dominant when ζ is increased. However, the large scale variations are the primary determinant of the aperture varia- tions for any positive ζ. A self-affine aperture field is analogous to fractional Brownian motion, if one views the spatial displacement |r | as a time-like variable and the change of aperture with displacement as equivalent to the spatial step of the Brownian motion. ζ< 1/2 With this analogy, corresponds to a subdiffusive or antipersistent process, ζ = 1/2 corresponds to a normal diffusive process, and 1/2 <ζ < 1 corresponds to a super diffusive process with ζ = 1 representing ballistic motion (Mandelbrot and Ness 1968). Interested readers are advised to consult the following references for more details on the theory of self-affine behavior and other applications, (Feder 1988; Peitgen and Saupe 1988; Falconer 2003). The algorithm used to generate spatially varying aperture fields, which was adapted from Méheust and Schmittbuhl (2001), uses random seeding and Fast Fourier Trans- forms to obtain isotropic self-affine fracture aperture fields on rectangular fractures. In order to generate circular fractures with a self-affine aperture field, a square fracture of length 2R was first generated and only the area of an inscribed circle was used for the analysis. Any shaped fracture could have been generated by first generating a large enough rectangular fracture to inscribe the desired shape of the fracture. Only the inscribed area formed by the arbitrary shaped fracture is used for the analysis. To obtain the desired σ of the fracture, the aperture value h is rescaled. h Fox et al. Geotherm Energy (2015) 3:21 Page 6 of 29 Real reservoirs are not static; there will be dynamic changes to the aperture and per- meability fields due to thermal, hydraulic, mechanical, and chemical (THMC) effects. In earlier works, these effects have been incorporated to model a reservoir’s thermal per - formance for both porous media (McDermott et al. 2006; Taron and Elsworth 2009) and discrete fractures (Ghassemi and Zhou 2011; Rawal and Ghassemi 2014). To achieve via- ble treatments, researchers have to specify quantitative deterministic models of THMC effects which themselves require further additional assumptions. The objective of our study was different and solely directed towards understanding how stochastic spatial aperture variations affect reservoir performance—not to incorporate coupled THMC effects. Ghassemi and Zhou (2011) modeled thermo-poroelastic stresses in a discrete fracture; their work showed that the area near the injection well experienced the largest increase in aperture, driven by local cooling that leads to a large temperature contrast compared to the initial reservoir temperature. In their modeling, the aperture near the injec- tion well increased from 0.23 to 0.38 mm over the operational lifetime of the reservoir. Regions in the middle of the reservoir and near the production well experienced less of an increase. Later, Rawal and Ghassemi (2014) considered silica dissolution, precipita- tion, and diffusion into the surrounding matrix along with thermo-poroelastic effects in a discrete fracture. They demonstrated that these effects changed the original constant aperture of the reservoir, with the cooler region near the injection well increasing by 10 % and the region near the production well decreasing by 24 %. The aperture changes in both of these studies are spatially smooth, which has the effect of lessening any strong perturbations to the base flow and the resulting temperature field. On the other hand, the randomly generated aperture fields with a self-affine description lack smoothness and areas of high and low aperture can occur anywhere in the reservoir, like the sam- ple apertures shown in Fig. 1. Spatial aperture fluctuations modeled as self-affine should not be considered a secondary effect with respect to THMC effects. The magnitude of the initial aperture variation is comparable with the magnitude of the changes due to THMC effects. While THMC induced aperture variations lead to systematic variations in pressure gradient and flow rate, the pre-existing aperture variations produce flow channeling. Fluid flow The fluid flow within the fracture aperture can be modeled as a flow of a Newtonian fluid in a thin channel. For fully developed laminar flow, if the aperture varies gradually, the Navier–Stokes equation is reduced to µ∇ v = ∇P − ρ g, (6) where P is the pressure, v is the velocity, ρ is the density of water, and µ is the dynamic viscosity. The conservation of mass can be expressed as � � ∇ · (ρ w v ) = m ˙ δ(r − r ). w l l (7) The summation on the right hand side refers to sources/sinks of mass, w is the local aperture or gap of the channel, and v is the fluid velocity averaged across the gap. The Fox et al. Geotherm Energy (2015) 3:21 Page 7 of 29 (a) (b) (c) (d) (e) Fig. 1 Dimensionless fracture surface temperature for the base case (a) and fractures exhibiting poor (b) and good (d) thermal performance (left hand side) with their respective aperture fields (c and e; right hand side). The aperture field is normalized by the mean aperture value w and the dimensionless aperture standard deviation σ /w and the Hurst coefficient ζ are equal to 0.25 and 0.8, respectively. The wellbore separation value is L/(2R) = 0.4 and the injection and production wells are located at (−0.4,0) and (0.4,0), respectively (marked by a black dot). All figures represent the same production period of 203 days, the time it takes for the base case to reach a dimensionless production temperature of  = 0.8 p,0 sources and sinks of mass in the fracture come from the connection to the wellbores, the injection (source) and production (sink) wells. Poiseuille flow describes pressure driven laminar flow in a channel with a no slip condition applied at the walls of the channel. For Poiseuille flow, the velocity averaged across the gap is �v� =− ∇P − ρ g , w (8) 12µ Fox et al. Geotherm Energy (2015) 3:21 Page 8 of 29 and this velocity is a function of the position r in the two-dimensional plane of the frac- ture. Combining Eqs. (7) and (8) yields ρ w −∇ · [∇P − ρ g] = m ˙ δ(r − r ), w l l (9) 12µ where the gradient operator only applies in the plane of the fracture. Often, Eq. (9) is referred to as the Reynolds equation and is derived using an approxi- mate lubrication analysis. As such, Eq. (9) is valid when Re(w/L ) ≪ 1 and w/L ≪ 1 x,y x,y (Deen 1998). L refers to the characteristic length scale of aperture variations in the x x,y and y direction and may be considered a “wavelength” of aperture variations. Here, the Reynolds number is defined as Re ≡ (ρ U w)/µ , where U is the characteristic velocity. w c c For a mass flow of 40 kg/s, a wellbore spacing L of 400 m, a width of 1 mm, and a char - acteristic velocity of U ∼˙m/(ρ wL), then Re ∼ 100, well below the critical Reynolds c w number for turbulence transition. The lubrication criteria Re(w/L ) ≪ and w/L ≪ 1 x,y x,y will also be satisfied when Re ∼ 100 because L > 1 m because large scale aperture var- x,y iations dominate when the aperture fluctuations are self-affine. Eq.  (9) has been exten - sively used to model fluid flow in fractures (Brown 1987; Zimmerman and Bodvarsson 1996; Méheust and Schmittbuhl 2001; Watanabe et  al. 2008). Because the aperture is raised to the third power, flow modeled using Eq.  (9) is sometimes said to follow the “cubic law”. The pressure field is first solved using Eq.  (9) with the finite element method (FEM). After solving for the pressure, the velocity is obtained from Eq.  (8). The injection and production wells are modeled as circles of radius 0.15  m to avoid singularities in pres- sure. The mass flux was evenly distributed along the circumference of the injection and production well. We specify a no mass flux condition on the boundary of the fracture and in the direction normal to the surface in order to model the fracture having a lim- ited areal extent with no leak off, that is the surrounding rock matrix is assumed to be impermeable. Heat transport The numerical model we have developed is a hybrid finite element and boundary ele - ment method for discretely fractured geothermal reservoirs. The three-dimensional heat conduction in the rock matrix is captured through a boundary element treatment while the fluid flow and thermal convection are solved via the FEM. As stated earlier, the rock matrix surrounding the discrete fractures is assumed to be impermeable with respect to fluid flow. The numerical model relies on 2D discretization of only the fracture surface and not the entire domain. The 2D discretization offers faster numerical run times over a 3D discretization of the entire reservoir and allows one to invest more computational effort on fracture scale details. The approach to use hybrid FEM and boundary integral representations is not new and our developed model was influenced by such works as Zhou et al. (2009), which studied the effect of thermoelastic stresses in a geothermal res - ervoir. Other hybrid numerical and analytical methods have been derived before, Cheng et al. (2001) used their developed methodology to model multi-dimensional heat effects in geothermal reservoirs and McDermott et al. (2009) used their hybrid model to solve Fox et al. Geotherm Energy (2015) 3:21 Page 9 of 29 for the flow, mass, and heat transport in a fracture, taking advantage of not having to dis - cretize a three dimensional domain. The fluid energy balance inside the fracture is ∂T ρ c w + ρ c w �v� · ∇T − k w∇ T = q . w w w w w s (10) ∂t The term q refers to the surface heat flux and c and k are the heat capacity and ther- s w w mal conductivity of water, respectively. By neglecting the transient heat storage in the liq- uid contained by the fracture and thermal conduction in the liquid (Pe ≡ U L/α ≫ 1 ), c w Eq. (10) becomes ρ c w�v� · ∇T = q . w w s (11) Unfortunately, the preceding equation cannot be solved independently because the surface heat flux is unknown and coupled to transient conduction of the rock mass sur - rounding the fracture. The surface heat flux is also related to the temperature of the sys - tem and can be represented by the 3D Green’s function and the intensity of the source. Here, the source is the fracture surface heat flux. For a source located at z and confined by the extent of the fracture, the temperature of the system is expressed as: t ′ ′ ′ q (x , y , t ) ′ ′ ′ ′ ′ ′ T (x, y, z, t) = T + G(x − x , y − y , z, t − t ) dx dy dt , r (12) ρ c r r 0 A ′ ′ q¯ (x , y , s) ′ ′ ′ ′ ¯ ¯ ¯ T (x, y, z, s) = T + G(x − x , y − y , z, s) dx dy , (13) ρ c r r where G is the 3D Green’s function for transient conduction and the overbar notation refers to the Laplace transform with s representing the Laplace space variable. The sub - script r refers to the thermophysical properties of the rock while T is the original res- ervoir rock temperature. The Laplace transform of the fracture heat balance equation is simply � � ρ c w v · ∇T =¯q . (14) w w s By working in the Laplace space, time discretization is explicitly avoided along with singularities at t = 0 and accuracy issues associated with early time discretization. Equa- tions (13) and (14) are solved simultaneously to obtain the surface heat flux and the tem - perature profile of the fracture at different time steps of interest. The inverse Laplace transform is computed using the Gaver–Stehfest algorithm (Stehfest 1970a, b) which determines the value of f(t) by computing several terms in a series involving f (s) at dif- ferent values of s. ′ ′ The matrix G (x − x , y − y , z) is computationally intensive to construct because every ij entry must be evaluated. The computational cost is O(N ) where N is the number of nodes. Computations can be greatly reduced by noticing that many of the entries in G ij ¯ ¯ are virtually zero and can be skipped in populating G . For determining G for a given ij ij s, entries that correspond to evaluating points that are farther than 20 α/s are skipped where α refers to the thermal diffusivity of the rock. The computational intensity can then be reduced to O(N log N ). The value 20 is arbitrary and can be changed to modify Fox et al. Geotherm Energy (2015) 3:21 Page 10 of 29 accuracy. This method is akin to a single term Fast Multipole Method as described by Greengard and Rokhlin (1997). For a 2500 node mesh, the computational time was reduced by a factor of about 20 with only a difference in the fifth significant figure for the outlet temperature. Alternatively, since the conduction is dominated by the direction perpendicular to the fracture, the 1D Green’s function can be used instead of the full 3D form. Doing so results in a less than 1 % error and over a 20-fold decrease in simulation time compared to using the truncation approximation of G . Additionally, the 1D simplification corre - ij sponds to a 400-fold decrease in simulation time compared to the fully populated G ij matrix. Fox et al. (2013) showed that 1D conduction is a reasonable approximation when β ≪ 1 β . Here, is a dimensionless parameter describing the ratio of the extent of ther- mal depletion normal to the plane of the fracture to the extent of depletion in the plane of the fracture. For dipole flow, β can be defined as β = 2k L/(mc ˙ ), where L is the spac- r w ing between the injection and production well. For typical flow conditions of 40 kg/s and 2 −4 a wellbore separation of 400 m, β = 2 × 10 . Nonetheless, for completeness, we used the full 3D heat conduction in this study. The dimensionless temperature in the reservoir is defined as � ≡ (T(x, y) − T ) /(T − T ), where T is the temperature of the fracture surface at point (x, y), T is the r w w temperature of the injected water and T is the initial rock temperature which remains constant in the far-field. The value of  ranges from zero to one. A value of zero refers to complete thermal exhaustion of the reservoir as the temperature in the fracture is now equal to the injected water temperature. If  = 1, there is no drawdown since the temperature in the fracture is still at the initial reservoir rock temperature. The variable refers to the dimensionless production temperature. If β ≪ 1 (1D conduction), then the dimensionless temperature is a function of the dimensionless wellbore spacing and dimensionless time: L t � = f , , (15) 2R τ 4 2 where τ ≡ L k ρ c /(mc ˙ ) . The dimensionless wellbore spacing is constrained to range r r r w from zero to one. The results for one mass flow rate and a set of times for a given dimen - sionless wellbore spacing are the same as for a different mass flow rate and a set of extrac - tion times corresponding to the same dimensionless times. By adding self-affine aperture fluctuations, the statistical values of the dimensionless temperature become functions of two additional parameters, σ /w and ζ. For this study, we have used a radius R of 500 m and an injected mass flow m ˙ of 40 kg/s. The thermophysical property values used were k = 2.4 W/(m °C), ρ = 2300  kg/m , c = 1000 J/(kg °C), and c = 4184 J/(kg °C). The r r r w value of τ is 2.1 × 10  s or 0.16 years for a wellbore spacing of 400 m and for the previ- ously listed values of the thermophysical properties and mass flow rate. Results and discussion Results for one fracture For an idealized single circular fracture with self-affine aperture variations, our mode - ling provides the statistical distribution of the dimensionless production temperatures as a function of wellbore spacing, duration of extraction, the variance of fracture aperture Fox et al. Geotherm Energy (2015) 3:21 Page 11 of 29 variations, and the Hurst coefficient. The extraction time for geothermal production selected for our study corresponded to a particular level of the thermal performance. In our study, this occurred when the dimensionless production temperature reached 0.8 for the case without aperture variations, which we call the base case. The variable  refers p,0 to the base case production temperature. The drawdown criteria will vary depending on the end-use of the thermal energy; direct thermal use operations can afford a greater drawdown compared to electricity production. However, using  = 0.8 is a realistic p,0 value for the drawdown criteria in many practical situations. For σ /w, ζ, and L/(2R), we chose a range of possible values for these parameters to model their impacts on thermal performance. The range of σ /w was chosen to be σ /w from 0.015 to 0.35. Granite core samples have been shown to have a value of of −4 0.42 (Hakami and Larsson 1996) but also low values around 5 × 10 (Sausse 2002). Because of the uncertainty in the standard deviation of the local aperture values and the difficulty in concluding how representative core sample values will be of practical in situ conditions in geothermal reservoirs, we believe the range chosen is adequate for simulating a large class of fractures. Above σ /�w� = 0.35, it becomes difficult to gener - ate fractures that do not have zero or negative aperture values. u Th s, the aperture reali - zations generated represent fractures that have been sufficiently propped open due to the combination of “self-propping” due to shear displacement and inflation as a result of higher fluid pressures within the fracture. As previously reported, the Hurst coeffi - cient has been observed to be always above 0.5 and is often reported to be around 0.8. However, the Hurst coefficient is a model specific parameter and is not often reported. It is plausible to think that there are instances where the Hurst coefficient is less than 0.5. Given its theoretical range of 0 to 1, the range of the Hurst coefficient chosen was from 0.1 to 1. The dimensionless wellbore spacing L/(2R) can vary anywhere from 0 to 1, depending on where the injection and production wells have been placed relative to the fracture size. In this study, we investigated the range from 0.2 to 0.8. We first present the results where the standard deviation of the aperture is varied. For the range of values investigated, five thousand fracture realizations were generated for each value of σ /w with L/(2R) and ζ kept constant at 0.4 and 0.8, respectively. In a second set of simulations, we repeated the process but now varying the Hurst coeffi - cient and keeping L/(2R) = 0.4 and σ /�w� = 0.25. Finally, in a third set, the wellbore spacing was varied and σ /w and ζ are set to 0.25, and 0.8, respectively. The statistical values of the thermal performance of most interest are the mean dimensionless produc- tion temperature  of the ensemble, the standard deviation of production values σ , p � and the fraction of cases with a higher production temperature than the base case F. Five-thousand was chosen to have an adequate sample size but also to be able to gener- ate results in a timely manner. Using ten-thousand fracture realizations did not result in any significant change in the derived results. Varying strength of aperture fluctuations σ /w For a range of from 0.015 to 0.35, five thousand fracture realizations were gener - ated for each value of σ /w. For each five thousand fracture realization data set, the values of  , σ , and F are reported based on the thermal performance for all realiza- p � tions at the time when the dimensionless production temperature reached 0.8 in the Fox et al. Geotherm Energy (2015) 3:21 Page 12 of 29 absence of aperture variations. A depiction of the fracture surface temperature for the base case, thermally under performing, and thermally over performing (compared with the base case) and their corresponding aperture fields is shown in Fig.  1. The dimension - less heat transfer area accessed at a certain time is defined as (1 − �) dA/A = 1 − ��� . u Th s, the mean fracture surface temperature  is a measure of the thermally accessed area of the fracture at a given extraction time. Comparing  for different fracture aper - ture distributions, a fracture with a lower value of  was able to access more area of the reservoir for heat transfer. For dipole flow, the majority of the flow occurs in the region surrounding the injection and production wells with less flow in the periphery of the fracture. The base case mean value for  is 0.816 when  = 0.8. For a fracture with p,0 smaller apertures in the region between the wells, the fluid is pushed to higher aper - ture regions in the periphery and the effective heat transport area increases, �� = 0.787 when the dimensionless production temperature reaches 0.8. As expected, this less channeled, beneficial flow field leads to a higher production temperature and slows the rate of thermal drawdown. On the other hand, when there exists a high aperture region in between the wells, the flow is further channeled from the injection to production well and the effective heat transfer area becomes narrower than the base case; the mean dimensionless fracture surface temperature is 0.854 when the dimensionless production temperature is 0.8. Consequently, the production well cools down rapidly. How the low and high aperture regions are arranged in the fracture is clearly important in determin- ing whether the reservoir will exhibit poor or good thermal performance. Figure 2 shows the distribution of thermal performance when σ /�w� = 0.25. Although a Gaussian curve does not fit the histogram perfectly, it does an adequate job of describ - ing the behavior. A Gaussian distribution of  is expected for a Gaussian aperture fluc - � � tuations with small amplitude variations. Inserting w + h into the governing equations for fluid velocity and fracture surface temperature and considering a regular perturba - 2 3 tion expansion for small ǫ = σ /�w�, one obtains � = � + ǫ� + ǫ � + O(ǫ ) for w, p 0 1 2 where  is the production temperature for the constant aperture case and  is a linear functional of the field h(r)/σ . A linear functional of a Gaussian field is a Gaussian vari - able. Figure 2 also reveals that we sometimes encounter enhanced thermal performance but more often encounter poor performing fractures with a dimensionless production temperature less than the base case. Figures  3 and 4 show how the strength of aper- ture variations affects the mean and standard deviation of the production temperature. As σ /w increases, the mean production temperature decreases quadratically and the standard deviation increases linearly. The trends reveal that as the aperture fluctuations increase, the thermal performance declines on average and is more variable. In effect, Fig.  4 says that predicting the reservoir’s thermal performance becomes more difficult with increasing values of σ / w . The linear behavior of Fig.  4 is expected due to the before mentioned linear behav- ior of  to small variations of h. In particular,  has zero standard deviation, while p 0 has an O(1) standard deviation so that the standard deviation of � ≈ � + ǫ� is 1 p 0 1 proportional to σ . The quadratic relationship between the mean production tempera - ture and the aperture standard deviation is also expected owing to the fact that  must ¯ ¯ ¯ be zero based on linearity. Thus, � − � = (σ /�w�) � , where  is a constant of 0 h 2 2 proportionality. Fox et al. Geotherm Energy (2015) 3:21 Page 13 of 29 Fig. 2 Histogram of the dimensionless production temperature  when the dimensionless wellbore spac- ing L/(2R) is equal 0.4, the dimensionless standard deviation of the aperture σ /w is equal to 0.25, and the Hurst coefficient ζ is equal to 0.8. The results are for the production period when the base case (no aperture variations) reaches  = 0.8. A Gaussian fit results in a mean and standard deviation of 0.783 and 0.072, p,0 respectively Fig. 3 Deviation of the mean dimensionless production temperature from the base case value  −  as p,0 p a function of the dimensionless standard deviation of the aperture σ /w. The Hurst coefficient ζ is equal to 0.8, the wellbore spacing L/(2R) is equal to 0.4, and the base case production temperature  is equal to 0.8 p,0 (reference value of no aperture variations). A quadratic fit is also shown Geothermal operations are sensitive to temperature drawdown and a premature drop in production temperature would reduce the economic attractiveness of the project. Each fracture in the ensemble has a lifetime corresponding to the time when  reaches 0.8. When σ /�w� = 0.1, which resulted in σ = 0.031, the 25 and 75 % quantiles of the distribution of the ratio of reservoir life to the reservoir life of the base case are equal to 0.83 and 1.16, respectively. For σ /�w� = 0.35, corresponding to σ = 0.092, the results h � are more drastic as the 25 and 75  % quantiles become 0.47 and 1.28, respectively. On the other hand, aperture variations are less important for lower values of σ /w; when � � σ / w = 0.05, the 25 and 75 % quantiles of reservoir life are closer together and equal to 0.91 and 1.08, respectively. Fox et al. Geotherm Energy (2015) 3:21 Page 14 of 29 Fig. 4 Standard deviation of the dimensionless production temperature σ for the ensemble of fracture realizations as a function of dimensionless standard deviation of the aperture σ /w when ζ = 0.8 and L/(2R) = 0.4. The production end time of interest is when the dimensionless production temperature  of the base case (no aperture variations) reaches 0.8. A linear fit is also shown For a normally distributed production temperature, the probability density function is equal to 1 (� − � ) p p f (� ) = exp . p (16) 2σ 2πσ From Fig. 2, it is reasonable to model the production temperature using a normal (Gauss- ian) distribution. The fraction of fractures that exhibit enhanced thermal performance is equal to one minus the integral of Eq. (16) with respect to  form 0 to  . The result is p p,0 1 � − � p,0 p F = erfc √ . (17) 2σ As discussed earlier, the mean and standard deviation of the production temperature distribution have a quadratic and a linear behavior, respectively: � − � = C , (18) p,0 p 1 � � and σ = C . � 2 (19) �w� Using Eqs. (18) and (19), the fraction of fractures that exhibit enhanced thermal per- formance compared to the base case can be expressed in terms of the standard deviation of the aperture fluctuations, 1 Kσ F = erfc √ , (20) 2�w� Fox et al. Geotherm Energy (2015) 3:21 Page 15 of 29 where K is the ratio of C /C . Since the error function is linear for small values of 1 2 K σ /( 2�w�), Eq. (20) can be approximated as 1 K σ F ≈ − √ . (21) 2 �w� 2π Figure  5 displays the linear decrease of F as a function of σ /w. The linear behavior of F results because the argument of the error function is small for the range studied (K = 1.03). The theoretical prediction for F matches well with the simulation results and is consistent with assuming that the distribution of  is Gaussian. The theoretical and simulated curves only deviate about 0.015 at the end of the range investigated. The frac - tion of fractures exhibiting adverse flow channeling was not observed to be a function of dimensionless time as F only varied temporally by less than 1 %. Thus, choosing a differ - ent dimensionless temperature stopping criteria will not change the reported values of F. The results shown indicate that the magnitude of aperture variations has a large effect on the thermal performance of the reservoir. For a specified standard deviation, the pro - duction temperature can be modeled as a Gaussian distribution. The fraction of frac - ture whose thermal performance benefits from a self-affine aperture field monotonically decreases with respect to the strength of the aperture variations and is never higher than 0.5. It ranged from 0.49 to 0.34 for the range of σ /w investigated. Thus, there are more ways to arrange isotropic self-affine aperture fluctuations that lead to a flow field which produces adverse rather than beneficial thermal performance. Varying the Hurst coefficient To investigate the effect of the Hurst coefficient ζ on thermal performance, simulations were performed for ten sets of five thousand fractures. Each set had a different value of ζ in the range from 0.1 to 1 and σ /w was kept constant at 0.25. Figure  6 displays 1D fracture aperture profiles for various values of ζ. Increasing ζ yields aperture fields with strong large scale correlations. In other words, the profile for ζ = 1 is smoother than the Fig. 5 The fraction of fractures exhibiting thermal performance enhancement F as a function of dimen- sionless standard deviation of the aperture σ /w. The Hurst coefficient ζ is equal to 0.8 and the wellbore spacing L/(2R) is equal to 0.4. A best fit line (dashed red) was added as well as the theoretical prediction (solid black) using a Gaussian distribution and the fitted parameters for the relationship of σ and � p Fox et al. Geotherm Energy (2015) 3:21 Page 16 of 29 Fig. 6 The fracture aperture h profile for various values of the Hurst coefficient ζ with the dimensionless standard deviation of aperture fluctuations σ /w kept at 0.25. For values of ζ of 0.1, 0.5, 0.8, and 1, the pro- files were shifted upward by w, 2w, and 3w to prevent them from overlapping profile for ζ = 0.1, which is characterized by many “peaks” and “valleys”. Table  1 sum- marizes the statistical results of the ensembles for each value of ζ, with a 95 % confidence interval derived by resampling based on the bootstrapping method (Efron and Tibshi- rani 1994). Table 1 indicates that the Hurst coefficient has little effect on the thermal performance given the confidence intervals of the statistical quantities. Zimmerman and Bodvarsson (1996) and Méheust and Schmittbuhl (2001) have shown that changing the Hurst coef- ficient had little affect on fluid flow. Since heat transport is mathematically an integral process that is dependent on large length scale fluid flow variations and even the vari - ance of the aperture is governed by the large scale, it is not surprising that varying the Hurst coefficient has little impact on the distribution of thermal performance. Of the statistical values from the ensemble, only σ changes appreciably, ranging from 0.063 to 0.073 with increasing ζ. If the aperture field is more correlated (larger value of ζ), the Table 1 Thermal performance with varying the Hurst coefficient ¯ F ζ σ 0.1 0.786 ± 0.002 0.063 ± 0.001 0.41 ± 0.01 0.2 0.785 ± 0.002 0.065 ± 0.001 0.40 ± 0.01 0.3 0.784 ± 0.002 0.067 ± 0.001 0.40 ± 0.01 0.4 0.784 ± 0.002 0.069 ± 0.001 0.40 ± 0.01 0.5 0.784 ± 0.002 0.069 ± 0.001 0.39 ± 0.01 0.6 0.784 ± 0.002 0.071 ± 0.001 0.40 ± 0.01 0.7 0.783 ± 0.002 0.072 ± 0.001 0.39 ± 0.01 0.8 0.783 ± 0.002 0.072 ± 0.001 0.39 ± 0.01 0.9 0.784 ± 0.002 0.073 ± 0.001 0.40 ± 0.01 1 0.784 ± 0.002 0.073 ± 0.001 0.40 ± 0.01 The statistical values are reported with a 95 % confidence interval determined using resampling based on the bootstrapping method Fox et al. Geotherm Energy (2015) 3:21 Page 17 of 29 region between the wells is more likely to have either a single peak or a single valley (low or high aperture regions). Having either a single peak or valley results in a more drastic difference in thermal performance than would be induced by the many peaks and valleys occurring for smaller ζ. The standard deviation of the production temperature does not change much for ζ> 0.6 because for these values the diameters of the peaks and valleys are likely to be equal or larger than the wellbore spacing. Varying wellbore separation The effect of self-affine aperture variations on different base flows can be understood by changing the value of the wellbore separation value, L/(2R). As L/(2R) gets smaller the streamlines becomes more varied and diverse. On the other hand, the streamlines become more uniform as the injection and production wells are pulled apart. For values of L/(2R) of 0.2, 0.4, 0.5, 0.6, and 0.8, the thermal performance of the five thousand frac - ture simulations for each value of L/(2R) were all performed for an end time when the respective L/(2R) base case production temperature reaches 0.8. To reach the  = 0.8 , the dimensionless time t/τ for L/(2R) of 0.2, 0.4, 0.5, 0.6, and 0.8 is equal to 4.7, 3.5, 2.8, 2.2, and 1.2 respectively. Although increasing the wellbore spacing increases the time to reach the drawdown criteria, the dimensionless time t/τ to reach the drawdown criteria drops for increasing values of wellbore spacing because τ is proportional to L . The val - ues of ζ and σ /w were kept constant at 0.8 and 0.25, respectively. Boxplots for the distribution of thermal performance for the different values of L/(2R) are plotted in Fig.  7. The boxplots show that the distribution of dimensionless produc - tion temperatures becomes more varied with increasing values of the wellbore spacing but peaks at L/(2R) = 0.4; the standard deviation of the dimensionless production tem- perature is greatest for this value. The mean value of the distribution and F decreases as L/(2R) decreases. Table 2 summarizes the statistical values of the data sets, with a 95 % confidence interval derived by resampling based on the bootstrapping method (Efron and Tibshirani 1994). The mean dimensionless production temperature is closest to the base case value when the wellbore spacing is small and there is less opportunity for the conductance to vary in the inter-wellbore region. The thermal performance distribution is most varied when L/(2R) = 0.4 because it corresponds to the typical length scale of one “peak” or “valley”, as can be see in Fig. 6. When the dimensionless wellbore spacing is smaller than 0.4, then the region in between the wells will have less aperture variations which results in a lower σ . The trend of lower values of F for larger values of wellbore spacing can be rationalized by considering the heat transfer area. The mean fracture surface tempera - ture  is a measure of the accessed area for heat transfer at a point in time; smaller val- ues of  mean more area has been accessed. The base case values for  are 0.95, 0.82, 0.74, 0.67 and 0.57 for L/(2R) = 0.2, 0.4, 0.5, 0.6, and 0.8, respectively. For larger wellbore spacing values, the base case is already harvesting energy from a large portion of the fracture and any deflection of the base flow caused by aperture variations will tend to decrease the heat transfer area as seen in the results for F, the fraction of fractures with better performance than the base case. In contrast, smaller wellbore spacing values, the base case heat transfer area is small and perturbations of the flow field are more likely to increase the heat transfer area. Fox et al. Geotherm Energy (2015) 3:21 Page 18 of 29 Fig. 7 Boxplot of dimensionless production temperatures for various values of wellbore spacing L/(2R). The dimensionless standard deviation σ /w is equal to 0.25 and the Hurst coefficient ζ is equal to 0.8. The lower and upper edge of the box represents the first and third quartiles while the solid red line represents the median value. The lower and upper whiskers of the boxplot are defined as 1.5 times the first and third quartile value, respectively, and the red plus signs represent the outliers. The dashed red line represents the base case (no aperture variations) value of 0.8 Table 2 Thermal performance with varying wellbore spacing L/(2R) L/(2R) σ F 0.2 0.793 ± 0.002 0.057 ± 0.001 0.42 ± 0.01 0.4 0.783 ± 0.002 0.072 ± 0.001 0.39 ± 0.01 0.5 0.779 ± 0.002 0.073 ± 0.001 0.38 ± 0.01 0.6 0.774 ± 0.002 0.068 ± 0.001 0.36 ± 0.01 0.8 0.763 ± 0.002 0.053 ± 0.001 0.28 ± 0.01 The statistical values are reported with a 95 % confidence interval determined using resampling based on the bootstrapping method Results for multiple fractures The previous section described the effect of a self-affine aperture field on the thermal performance when there is only one fracture. A relevant situation to further model would be a reservoir consisting of two fractures. We consider a reservoir consisting of two fractures with equal radii that do not intersect with one another, but rather are hydraulically connected via the injection and production well. The pressure drop along the wellbore is assumed negligible because w/r ≪ 1, where the wellbore radius r is around 10 cm and the aperture is around 1 mm. The mass flow was doubled to 80 kg/s, but the value of m ˙ /A remains the same as that in the one fracture system because the total surface area of the system is twice as large. With two fractures the flow field exhibits both inter- and intrafracture flow variations. As has been seen in the one fracture reservoir, local aperture variations disturb the base- flow creating intrafracture flow variations. However, aperture variations can also modify the ratio of the mean flow rate between fractures through altering the flow impedance or resistance to mean flow of each fracture and affect the reservoir’s overall pressure drop. The flow impedance of the fractures in the reservoir is what drives the global distribu - tion of flow in the reservoir. Fox et al. Geotherm Energy (2015) 3:21 Page 19 of 29 Since the pressure drop P is linear with the flow rate, P ∝˙m, the flow split fractions for fracture A and B are φ = , A (22) I + I A B φ = 1 − φ , B A (23) where I ≡ �P/m ˙ and represents the flow impedance of fracture i. The system can be i i viewed as two “resistors” in parallel where the pressure drop is the same between the two fractures. The above equation demonstrates how fracture aperture variations mod - ify the interfracture flow. Since it does not matter whether a fracture is labeled fracture A or B, the labels used in Eqs. (22) and (23) are arbitrary and we may define the variable φ as max(φ , φ ), which ranges from 0.5 to 1. A B The hydraulic aperture is another way to view the varying pressure drop for each frac - ture. In terms of the pressure drop of a fracture P, the hydraulic aperture is defined as 1/3 �P w ≡ w , (24) H 0 �P where P is the reference pressure drop in a smooth fracture for an aperture value of w receiving the same mass flux. The cube root is a result of the cubic law described in Eq. (9). The hydraulic aperture is the aperture value needed in a smooth fracture to have the same pressure drop in a fracture with aperture variations. The hydraulic aperture is a useful metric in quantifying how spatial aperture variations control the pressure drop and affect interfracture flow variations. The effects of fracture aperture variations in two and three fracture systems was stud - ied by determining the thermal performance of an ensemble of fracture realizations. The wellbore spacing and the Hurst coefficient were kept constant and the standard devia - tion of the aperture was varied over the same range considered before. Finally, we exam- ined the efficacy of flow control methods in minimizing any negative effects aperture variations have on the thermal performance of multifracture systems. Varying strength of aperture fluctuations Two thousand five hundred fracture pairs with L/(2R) = 0.4 and ζ = 0.8 were generated and their thermal performance was simulated. Reassigning the pairing had an insignifi - cant effect on the results presented. Figures  8 and 9 show the histogram of the produc- tion temperature and the flow split fraction when σ /�w� = 0.25. It is apparent that the flow split can be modeled as a Gaussian distribution and compared to Fig.  2, the thermal performance is closer to Gaussian than that of the one fracture system. Figure  10 shows that larger φ values, corresponding to more unequal flow splits, are likely to result in lower production temperatures. Also, plotted in the figure is the pro - duction temperature as a function of the flow split if the fractures had no surface vari - ations. In the absence of aperture variations,  is always monotonically decreasing with respect to φ. There are simulation outcomes that lie above the smooth case curve, but systems that exhibit a boost in thermal performance for larger flow split fractions Fox et al. Geotherm Energy (2015) 3:21 Page 20 of 29 Fig. 8 Histogram of the thermal performance of the two fracture system when L/(2R) = 0.4, σ /�w� = 0.25, and ζ = 0.8. A Gaussian fit results in a mean and standard deviation of 0.751 and 0.0625, respectively � � Fig. 9 Histogram of the fractional flow split of the two fracture system when L/(2R) = 0.4, σ / w = 0.25, and ζ = 0.8. By definition the mean is 0.5 and standard deviation is 0.156, respectively become increasingly rare. In general, having a flow split deviate from equality is likely to result in undesirable thermal performance. As plotted in Fig. 11, the uncertainty in ther- mal performance has a linear relationship with the uncertainty in φ. As with the single fracture reservoir, the fraction of reservoirs that benefit from aperture variations decreases with increasing mean square of aperture variations. The fraction of reservoirs with a production temperature higher than the base case, F, as a function of the dimensionless standard deviation, is plotted in Fig.  12. For comparison, the single fracture results are plotted in the same figure. The theoretical fit was obtained by using Eq. (20) along with the fits of and . Compared to the single fracture case, � p the values of both σ and  are smaller. The curve is no longer linear since K is now � p equal to 3.00 and as a result, the linear approximation of the error function breaks down. Fox et al. Geotherm Energy (2015) 3:21 Page 21 of 29 Fig. 10 Dimensionless production temperature as a function of fractional flow split when L/(2R) = 0.4, � � σ / w = 0.25, and ζ = 0.8. The solid black line represents the production temperature for smooth fractures when the flow split varies. One and two standard deviations from the mean are plotted to help visualize how the data is distributed Fig. 11 The standard deviation of the dimensionless production temperature σ as a function of the stand- ard deviation of the fraction flow split σ when L/(2R) = 0.4 and ζ = 0.8. A linear fit is also shown A portion of the decreased thermal performance in the two fracture reservoir results solely from the existence of an unequal flow split between the fractures. The reduction in the production temperature for two smooth fractures with increasing flow split ratio is indicated by the solid line in Fig. 10 and we will refer to this function as g(φ). A metric of the additional loss of thermal performance in the two fracture reservoir compared with the case of two smooth fractures with unequal apertures is F , which we define as the fraction of two fracture reservoirs with � > g(φ), i.e., the fraction that beat the smooth fracture reservoir with the same flow split. The results for F are higher than F because the smooth case thermal performance drops with respect to increasing φ. How- ever, F is closer to the results of F for the two fracture results than to the single fracture results. These results demonstrate that the primary determinant of the diminished per - formance of the two fracture system is the aperture variation and flow channeling within Fox et al. Geotherm Energy (2015) 3:21 Page 22 of 29 Fig. 12 The fraction of two fracture system exhibiting thermal performance enhancement as a function of σ /w. The Hurst coefficient ζ is equal to 0.8 and the wellbore spacing L/(2R) is equal to 0.4. The curves give theoretical predictions based on a Gaussian distribution of  . The squares are the fraction of reservoirs that outperform a base case with an equal flow split, F. The diamonds are the fraction of reservoirs that outper- form a base case with a constant aperture within each fracture but the same flow split as observed in the simulations, F . The one fracture results and the theoretical fit was also plotted each reservoir rather than being simply a result of the difference of hydraulic aperture of the two fractures. The reduced thermal performance compared to a one fracture system can be explained by the relationship between pressure drop and the thermal performance in an individ- ual fracture. In general, a fracture with a small aperture region in between the injec- tion and production wells will exhibit enhanced thermal performance because the flow will expand and sweep through the reservoir more effectively. The reverse is true; a high aperture region in between the wells will promote flow short circuiting and the flow will primarily harvest thermal energy in a narrow area. This behavior was illustrated in Fig.  1, where the temperature of the fracture surface was plotted along side their respective aperture fields. Furthermore, a low aperture wellbore region requires a higher pressure drop to push the fluid through the constricted region surrounding both the injection and production wells. The covariance of w and  for the single fracture is −0.65 when H p L/(2R) = 0.4, ζ = 0.8, and σ /�w� = 0.25, which supports our qualitative argument. For single fracture reservoirs, poor performing fractures exhibit a smaller pressure drop. However, in a two fracture reservoir, the poor thermally performing reservoir will receive a larger portion of the flow and the reservoir with good thermal performance receives a smaller portion of the flow. Thus, it is not simply the existence of an unequal flow split but the fact that the higher flow goes preferentially to the fracture with more channeling that accounts for the poorer performance of two fracture reservoirs com- pared with single fracture reservoirs. However, the negative impact on thermal perfor- mance can be mitigated if the flow split between the reservoirs can be manipulated. Three fracture system We can extend the analysis to a three fracture system, following the same method used for the two fracture reservoir. The system is treated as three resistors in parallel, Fox et al. Geotherm Energy (2015) 3:21 Page 23 of 29 all with the same pressure drop. Plotted in Fig.  13 are the fitted Gaussian probability distribution of the production temperature for one, two, and three fracture reservoir systems for when L/(2R) = 0.4, σ /�w� = 0.25, and ζ = 0.8. With the addition of more fractures, the distribution of thermal performance becomes less varied and the mean shifts to lower values. As expected, the thermal performance of the system gets worse for increasing values of the mean square of aperture fluctuations. For the range of σ /�w� = 0.015 − 0.35 , F monotonically decreases from 0.49 to 0.074. The three frac - ture system has the worst thermal performance among the cases investigated. Using the same argument as the two fracture system, it is more likely that a low impedance frac- ture has individually poor thermal performance and diverts most of the fluid into that fracture. With three fractures, the probability of having the aforementioned outcome increases and thermal performance suffers to a greater extent than the one or two frac - ture reservoirs. Flow control To mitigate the adverse impact of fracture surface variations on the interfracture flow variations, flow control methods can modify the flow split and reduce the negative effects associated with flow split disparity. In a real reservoir, flow control is achieved by using interflow control devices (ICD). An ICD works by increasing the resistance of flow from the wellbore to the formation or fracture. The increased resistance is achieved using nozzles, orifices, and helical channels on the base pipe. For more information on ICD, please refer to Al-Khelaiwi et al. (2010) and Birchenko et al. (2010). Figure  14 shows the boxplot of the dimensionless production temperature distribu- tion of the two fracture system when L/(2R) = 0.4, ζ = 0.8, and σ /�w� = 0.25 for three cases: no flow control (untreated), 50/50 flow split, and the optimal flow split. A 50/50 flow split was chosen as a good guess on how to best distribute the flow. The optimal Fig. 13 Fitted Gaussian probability density function of dimensionless production temperature for one, two, and three fracture system when L/(2R) = 0.4, σ /�w� = 0.25, ζ = 0.8, and with the base case value of = 0.8. The Gaussian distributions were chosen to match the mean and standard deviation of the simu- p,0 lated temperature results Fox et al. Geotherm Energy (2015) 3:21 Page 24 of 29 Fig. 14 Boxplot of dimensionless production temperatures for three cases: no flow control (untreated), 50/50 flow split, and the optimal flow split. The dimensionless standard deviation of the aperture σ /w is kept at 0.25, the Hurst coefficient ζ is equal to 0.8, and the wellbore spacing L/(2R) is equal to 0.4. The lower and upper edge of the box represents the first and third quartiles while the solid red line represents the median value. The lower and upper whiskers of the boxplot are defined as 1.5 times the first and third quartile value, respec- tively, and the red plus sign represents the outliers flow split represents the value of φ that optimizes the production temperature. In prac- tice, the optimal flow split would be very difficult to predict, as one would first need to know the structure of the reservoir from solving an inverse problem based on tracer tests. Nonetheless, it is calculated to be used as an upper bound on thermal perfor- mance. The boxplot of the 50/50 and the optimal flow scenario show that there is only a marginal increase in thermal performance when the optimal value is used. The 25, 50 % (median), and 75 % quantiles of the distribution for the change in production tempera- ture for the equal split are 0.0023, 0.0149, and 0.0458, respectively. For the optimal split, the 25, 50  % (median), and 75  % quantiles of the distribution are 0.0043, 0.0194, and 0.0554, respectively. Another way to view the flow control results is to plot the change in dimensionless production temperature  versus the flow split fraction. For the 50/50 flow split sce - nario, the results for  versus φ are plotted in Fig.  15. The equal flow split strategy does result in worse thermal performance 13 % of the time compared to no flow control. However, 75  % of these cases only saw a minimal drop in the production temperature, less than 0.004. Additionally, using the optimal flow split in these cases led to an increase in  of no more than 0.0019 with improvement occurring in only 75 % of the cases. The figure reveals that as the flow split fraction moves away from a 50/50 split, the improve - ments of applying flow control increases; the covariance of φ and  is 0.88. When the fracture surface is less varied, the distribution of thermal performance is also less varied. Consequently, flow control will become less important. When σ /�w� = 0.05, 75  % of the cases only had a  of less than 0.002 when the optimal value of φ was used. Compared to the F value for the one fracture reservoir, applying flow control for the range of σ /w investigated was capable of increasing F to a range of 0.50–0.31, which represents 100–89 % of the value of the one fracture reservoir. Plotted in Fig. 16 are the F values for the 50/50 flow split as well as the values for the one and two fracture cases (no flow control). It is only when σ /w > 0.25 that F of the equal flow split is appreciably lower than the one fracture reservoir. Fox et al. Geotherm Energy (2015) 3:21 Page 25 of 29 Fig. 15 The change of production temperature  , using an equal flow split strategy as a function of the original flow split fraction when L/(2R) = 0.4, σ /�w� = 0.25, and ζ = 0.8. The original flow split is normally distributed and one and two standard deviations from the mean are plotted to help visualize how  is distributed The fact that flow control was able to bring F near the single fracture results indicates that interfracture flow variations inhibit thermal performance more than intrafracture flow variations. The results in Fig.  12 showed that these detrimental interfracture flow variations consist of flow preferentially going to poor thermal performing fractures rather than simply involving unequal flow splits. Nonetheless a comparison of our two flow control strategies indicates that an effort to divert flow preferentially away from fractures exhibiting channeling yields only a marginal improvement over a flow control strategy that yields an equal flow split. The flow control results reveal that for equally sized fractures, a 50/50 flow split is good practice. Knowing the optimal split does not result in a large improvement in ther- mal performance. Flow control would be limited if the reservoir consists of fracture- to-fracture intersections as ICD can only control flow in the fractures that are directly Fig. 16 The fraction of two fracture system exhibiting thermal performance enhancement as a function of σ /w for the case of equal flow split. The Hurst coefficient ζ is equal to 0.8 and the wellbore spacing L/(2R) is equal to 0.4. The results for the one fracture, and untreated two fracture cases are also plotted Fox et al. Geotherm Energy (2015) 3:21 Page 26 of 29 connected to the wellbore. Fracture-to-fracture intersections, therefore, open up new possible flow paths rather than just partitioning the flow between two independent flow paths where the path with the least impedance takes the majority of the flow. Thus, the presence of fracture-to-fracture intersections may not result in a diminished thermal performance, as was shown in the case of multiple fractures only hydraulically con- nected via the wellbores. Conclusions Isotropic self-affine fracture aperture variations limit a reservoir’s thermal performance. There are more ways to arrange aperture variations that causes flow paths that lead to premature drawdown than flow pathways that enhance heat transfer. We have shown σ /w that the dimensionless standard deviation of the aperture is more important in controlling thermal performance than the Hurst coefficient ζ. Specifically, both the uncertainty in thermal performance, σ , and the number of fracture realizations that benefit from aperture variations have linear relationships to the root mean square of the aperture variations. However, the effects of aperture variations become nearly insignifi - cant when σ /w < 0.1. For the investigated range of σ /w from 0.015 to 0.35, 34 to h h 49  % of the fracture realizations benefited from aperture variations. Reservoirs with a larger wellbore separation have more uniform base flow, flow in the absence of aper - ture variations, that is more negatively affected by aperture variations. A reservoir with a large wellbore separation already accesses a larger area of the fracture for heat exchange and any flow channeling induced by aperture variations will more likely decrease the effective heat exchange area. Our finding that aperture variations will lead to adverse thermal performance, more often than not, is in agreement with the previous work of Neuville et al. (2010a, b), despite the fact that their study did not couple heat transport in the fracture with that in the surrounding rock and used a uniform base flow in a rectan - gular fracture. Additionally, our use of statistical values for an ensemble of realizations, including the mean and standard deviation of production temperatures, offers a better and more practical understanding of the effects fracture aperture variations have on dis - cretely fractured geothermal reservoirs. We also demonstrated that fracture aperture variations alter the interfracture flow split in a system with multiple parallel fractures. The number of reservoirs exhibiting diminished thermal performance increased compared a reservoir with one fracture and the same mass flow per fracture area. Simply having a higher flow rate in one fracture will more than likely lead to a worse thermal performance because the fracture receiving most of the flow will experience thermal breakthrough more quickly. Further deteriora - tion in thermal performance arises because fractures with a larger flow impedance usu - ally have better thermal performance but receive a smaller portion of the flow. However, applying flow control can alleviate the problems with undesirable flow splits, suggesting a need for inflow control devices specifically suited for geothermal applications rather than “off the shelf ” devices made for oil and natural gas operations. Our analysis suggests that a major cause of poor thermal performing reservoirs is due to interfracture flow variations. With more fractures, the reservoir is more susceptible to adverse thermal performance. In general, more fractures in the reservoir lead to a drop in the variability of thermal performance but also leads to a higher likelihood of having Fox et al. Geotherm Energy (2015) 3:21 Page 27 of 29 a worse production temperature compared to the base case scenario of no aperture variations. The multiple fracture systems simulated in this study were only hydraulically connected via the injection and production well. Real EGS reservoirs, however, like the one at Soultz-sous-Forêts will have fractures that intersect with each other (Sausse et al. 2010). With fracture-fracture intersections, flow control devices become less effective because they can only control well-to-fracture flow. The present study dealt only with a static reservoir where the spatially varying aper - ture field stayed constant with time. Real reservoirs, however, will experience ther - mal, hydraulic, mechanical, and chemical (THMC) changes that alter the aperture field throughout the energy extraction process. In addition to THMC effects, fracture shearing/slipping during heat extraction will also change the aperture field. In general, aperture change due to THMC effects will increase permeability/aperture in the region enclosed by the injection and production well (Taron and Elsworth 2009; Ghassemi and Zhou 2011). Fracture slippage will lead to an increase in the standard deviation of the aperture fluctuations and result in preferential flow channeling in the direction perpen - dicular to the shear displacement (Yeo et  al. 1998; Plouraboué et  al. 2000; Watanabe et al. 2008). Although our work did not incorporate fracture slippage and THMC effects it, nonetheless, provides a framework through which to understand how spatial aperture fluctuations affect thermal performance. It shows that fracture surface variations are a reservoir characteristic that researchers should consider in modeling discretely fractured geothermal reservoirs and performing techo-economic analysis (Beckers et  al. 2014; Held et al. 2014). Our systematic approach for generating spatially varying apertures can be used in conjunction with modeling the effects that THMC and fracture slippage have on reservoir thermal performance, by providing a base aperture field that already has a spatially varying aperture field. The developed numerical framework and the lessons learned from this study are cur - rently being applied to a field laboratory at the Altona Flat Rocks in Plattsburgh, New York. The Altona reservoir represents a meso-scale field site, larger than a bench scale, artificial reservoir but smaller than a full sized geothermal reservoir. The Altona reser - voir consists of Cambrian Potsdam Sandstone formation, where a five spot well system has been drilled in a 10 m by 10 m area, intersecting an isolated horizontal fracture 7.6 m deep (Becker and Tsoflias 2010). Ground penetrating radar has been used to determine the spatial aperture of the fracture (Tsoflias and Becker 2008), creating a reservoir with a characterized subsurface. The site has been used to conduct tracer tests, demonstrating the presence of flow channelizations caused by the spatially varying aperture field (Cast - agna et al. 2011; Becker et al. 2015; Hawkins et al. 2015). The site resembles the reservoir model used in this study, with an induced dipole flow, Gaussian distribution of aper - ture, and a self-affine spatial correlation. Future experimentation at the site will include deployment of thermally degrading and sorbing tracers and reservoir heating through injection of hot water as an analog of the geothermal energy extraction process. Authors’ contributions This study was performed by DBF with advice on the planning, design, and geothermal application of the model pro- vided by DLK and JWT. DBF wrote the manuscript and it was revised and the final version approved by all three authors. All authors read and approved the final manuscript. Author details 1 2 Chemical and Biomolecular Engineering, Cornell University, Ithaca 14850, USA. Cornell Energy Institute, Cornell University, Ithaca 14850, USA. Fox et al. Geotherm Energy (2015) 3:21 Page 28 of 29 Acknowledgements The authors would like to thank the National Science Foundation IGERT program (contract 0966045), the U.S. Depart- ment of Energy (contract DE-EE0000842), and the Cornell Energy Institute for partial financial support for this work. They would also like to thank their colleagues at Cornell University, specifically Koenraad Beckers, Adam Hawkins, Maciej Łukawski, and Calvin Whealton for their input and suggestions about this work. The authors appreciate the feedback and advice received from the two anonymous reviewers. Their input helped strengthen this work. Competing interests The authors declare that they have no competing interests. Received: 5 August 2015 Accepted: 30 September 2015 References Al-Khelaiwi FT, Birchenko VM, Konopczynski MR, Davies D. Advanced wells: a comprehensive approach to the selection between passive and active inflow-control completions. SPE Prod Oper. 2010;25:305–26. Auradou H, Drazer G, Boschan A, Hulin J-P, Koplik J. Flow channeling in a single fracture induced by shear displacement. Geothermics. 2006;35(5–6):576–88. Auradou H, Boschan A, Chertcoff R, DAngelo, M-V, Hulin J-P, Ippolito I. Miscible transfer of solute in different model frac- tures: From random to multiscale wall roughness. Comptes Rendus Geoscience. 2010;342(7–8):644–52. Becker MW, Tsoflias GP. Comparing flux-averaged and resident concentration in a fractured bedrock using ground pen- etrating radar. Water Resour Res. 2010;46(9):1–12 Becker MW, Tsoflias GP, Baker M, Hawkins A. Confirmation of hydraulic, tracer, and heat transfer characterization of a frac- tured bedrock using ground penetrating radar. In: Proceedings of the Fortieth Workshop on Geothermal Reservoir Engineering, Stanford University, January 26–28, 2015, Stanford, CA. 2015. Beckers KF, Lukawski MZ, Anderson BJ, Moore MC, Tester JW. Levelized costs of electricity and direct-use heat from Enhanced Geothermal Systems. J Renew Sustain Energy. 2014;6(1):1–15 Birchenko VM, Muradov KM, Davies DR. Reduction of the horizontal well’s heeltoe effect with inflow control devices. J Pet Sci Eng. 2010;75:244–50. Boffa JM, Allain C, Chertcoff R, Hulin JP, Plouraboué F, Roux S. Roughness of sandstone fracture surfaces: Profilometry and shadow length investigations. Eur Phys J B. 1999;7(2):179–82. Bouchaud E. Scaling properties of cracks. J Phys Condens Matter. 1997;9(21):4319–44. Brown SR. Fluid flow through rock joints: the effect of surface roughness. J Geophys Res Solid Earth. 1987;92(B2):1337–47. Brown SR, Scholz CH. Broad bandwidth study of the topography of natural rock surfaces. J Geoph Res Solid Earth. 1985;90(B14):12575–82. Brown SR, Kranz RL, Bonner BP. Correlation between the surfaces of natural rock joints. Geophy Res Lett. 1986;13(13):1430–3. Bruderer-Weng C, Cowie P, Bernabé Y, Main I. Relating flow channelling to tracer dispersion in heterogeneous networks. Adv Water Resour. 2004;27(8):843–55. Bruel D, Cacas MC, Ledoux E, de Marsily G. Modelling storage behaviour in a fractured rock mass. J Hydrol. 1994;162(3–4):267–78. Cacas MC, Ledoux E, de Marsily G, Tillie B, Barbreau A, Durand E, Feuga B, Peaudecerf P. Modeling fracture flow with a sto - chastic discrete fracture network: Calibration and validation 1. The flow model. Water Resour Res. 1990;26(3):479–89. Castagna M, Becker MW, Bellin A. Joint estimation of transmissivity and storativity in a bedrock fracture. Water Resour Res. 2011;47(9):1–19. Cheng AHD, Ghassemi A, Detournay E. Integral equation solution of heat extraction from a fracture in hot dry rock. Int J Numer Anal Methods Geomech. 2001;25(13):1327–38. Deen WM. Analysis of Transport Phenomena. New York: Oxford University Press; 1998. Drazer G, Koplik J. Permeability of self-affine rough fractures. Phys Rev E. 2000;62(6):8076–85. Drazer G, Koplik J. Transport in rough self-affine fractures. Phys Rev E. 2002;66:026303. Efron B, Tibshirani RJ. An Introduction to the Bootstrap. New York: Chapman & Hall; 1994. Falconer K. Fractal geometry: mathematical foundations and application. 2nd ed. Chichester: Wiley; 2003. Feder J. Fractals physics of solids and liquids. New York: Plenum Press; 1988. Fox DB, Sutter D, Beckers KF, Lukawski MZ, Koch DL, Anderson BJ, Tester JW. Sustainable heat farming: Modeling extrac- tion and recovery in discretely fractured geothermal reservoirs. Geothermics. 2013;46:42–54. Ghassemi A, Zhou X. A three-dimensional thermo-poroelastic model for fracture response to injection/extraction in enhanced geothermal systems. Geothermics. 2011;40(1):39–49. Glover PWJ, Matsuki K, Hikima R, Hayashi K. Fluid flow in synthetic rough fractures and application to the Hachimantai geothermal hot dry rock test site. J Geophys Res Solid Earth. 1998;103(B5):9621–35. Greengard L, Rokhlin V. A fast algorithm for particle simulations. J Comput Phys. 1997;135(2):280–92. Hakami E, Larsson E. Aperture measurements and flow experiments on a single natural fracture. Int J Rock Mech Min Sci Geomech Abstr. 1996;33(4):395–404. Hawkins A, Fox D, Zhao R, Tester J, Cathles L, Koch D, Becker M. Predicting thermal breakthrough from tracer tests: Simulations and observations in a low temperature field laboratory. In: Proceedings of the Fortieth Workshop on Geothermal Reservoir Engineering, Stanford University, January 26–28, 2015, Stanford, CA. 2015. Held S, Genter A, Kohl T, Koelbel T, Sausse J, Schoenball M. Economic evaluation of geothermal reservoir performance through modeling the complexity of the operating EGS in Soultz-sous-Forêts. Geothermics. 2014;51:270–80. Jupe AJ, Bruel D, Hicks T, Hopkirk R, Kappelmeyer O, Kohl T, Kolditz O, Rodrigues N, Smolka K, Willis-Richards J, Wallroth T, Xu S. Modelling of a European prototype HDR reservoir. Geothermics. 1995;24(3):403–19. Fox et al. Geotherm Energy (2015) 3:21 Page 29 of 29 Kolditz O. Modelling flow and heat transfer in fractured rocks: Dimensional effect of matrix heat diffusion. Geothermics. 1995;24(3):421–37. Mandelbrot BB, Ness JWV. Fractional Brownian motions, fractional noises and applications. SIAM Rev. 1968;10(4):422–37. McDermott CI, Randriamanjatosoa ARL, Tenzer H, Kolditz O. Simulation of heat extraction from crystalline rocks: the influ- ence of coupled processes on differential reservoir cooling. Geothermics. 2006;35(3):321–44. McDermott CI, Walsh R, Mettier R, Kosakowski G, Kolditz O. Hybrid analytical and finite element numerical modeling of mass and heat transport in fractured rocks with matrix diffusion. Comput Geosci. 2009;13(3):349–61. Méheust Y, Schmittbuhl J. Geometrical heterogeneities and permeability anisotropy of rough fractures. J Geophys Res Solid Earth. 2001;106(B2):2089–102. Méheust Y, Schmittbuhl J. Scale effects related to flow in rough fractures. Pure Appl Geophys. 2003;160(5–6):1023–1050. Neuville A, Toussaint R, Schmittbuhl J. Fracture roughness and thermal exchange: a case study at Soultz-sous-Forêts. Comptes Rendus Geosci. 2010a;342(7–8):616–25. Neuville A, Toussaint R, Schmittbuhl J. Hydrothermal coupling in a self-affine rough fracture. Phys Rev E. 2010b;82:036317. Nicol DAC, Robinson BA. Modelling the heat extraction from the Rosemanowes HDR reservoir. Geothermics. 1990;19(3):247–57. Nordqvist AW, Tsang YW, Tsang CF, Dverstorp B, Andersson J. A variable aperture fracture network model for flow and transport in fractured rocks. Water Resour Res. 1992;28(6):1703–13. Peitgen HO, Saupe D. The Science of Fractal Images. New York: Springer; 1988. Plouraboué F, Kurowski P, Hulin J-P, Roux S, Schmittbuhl J. Aperture of rough cracks. Phys Rev E. 1995;51(3):1675. Plouraboué F, Hulin J-P, Roux S, Koplik J. Numerical study of geometrical dispersion in self-affine rough fractures. Phys Rev E. 1998;58:3334–46. Plouraboué F, Kurowski P, Boffa J-M, Hulin J-P, Roux S. Experimental study of the transport properties of rough self-affine fractures. J Contam Hydrol. 2000;46(3–4):295–318. Ponson L, Auradou H, Pessel M, Lazarus V, Hulin J-P. Failure mechanisms and surface roughness statistics of fractured fontainebleau sandstone. Phys Rev E. 2007;76(3):036108. Rawal C, Ghassemi A. A reactive thermo-poroelastic analysis of water injection into an enhanced geothermal reservoir. Geothermics. 2014;50:10–23. Sausse J. Hydromechanical properties and alteration of natural fracture surfaces in the Soultz granite (Bas-Rhin, France). Tectonophysics. 2002;348(1–3):169–85. Sausse J, Dezayes C, Dorbath L, Genter A, Place J. 3D model of fracture zones at Soultz-sous-Forêts based on geo- logical data, image logs, induced microseismicity and vertical seismic profiles. Comptes Rendus Geosci. 2010;342(7–8):531–45. Schmittbuhl J, Gentier S, Roux S. Field measurements of the roughness of fault surfaces. Geophys Res Lett. 1993;20(8):639–41. Schmittbuhl J, Schmitt F, Scholz C. Scaling invariance of crack surfaces. J Geophys Res Solid Earth. 1995;100(B4):5953–73. Stehfest H. Algorithm 368: Numerical inversion of Laplace transforms [d5]. Commun ACM. 1970a;13(1):47–9. Stehfest H. Remark on algorithm 368: Numerical inversion of Laplace transforms. Commun ACM. 1970b;13(10):624. Taron J, Elsworth D. Thermal-hydrologic-mechanical-chemical processes in the evolution of engineered geothermal reservoirs. Int J Rock Mech Min Sci. 2009;46(5):855–64. Tsang YW, Tsang CF. Flow channeling in a single fracture as a two-dimensional strongly heterogeneous permeable medium. Water Resour Res. 1989;25(9):2076–80. Tsoflias GP, Becker MW. Ground-penetrating-radar response to fracture-fluid salinity: Why lower frequencies are favorable for resolving salinity changes. Geophysics. 2008;73(5):25–30. Watanabe N, Hirano N, Tsuchiya N. Determination of aperture structure and fluid flow in a rock fracture by high-resolution numeri- cal modeling on the basis of a flow-through experiment under confining pressure. Water Resour Res. 2008;44(6):1–11. Watanabe N, Wang W, McDermott CI, Taniguchi T, Kolditz O. Uncertainty analysis of thermo-hydro-mechanical coupled processes in heterogeneous porous media. Comput Mech. 2010;45(4):263–80. Yeo IW, de Freitas MH, Zimmerman RW. Eec ff t of shear displacement on the aperture and permeability of a rock fracture. Int J Rock Mech Min Sci. 1998;35(8):1051–70. Zhou XX, Ghassemi A, Cheng AH-D. A three-dimensional integral equation model for calculating poro- and thermoe- lastic stresses induced by cold water injection into a geothermal reservoir. Int J Numer Anal Methods Geomech. 2009;33(14):1613–40. Zimmerman RW, Bodvarsson GS. Hydraulic conductivity of rock fractures. Transp Porous Media. 1996;23(1):1–30.

Journal

Geothermal EnergySpringer Journals

Published: Oct 28, 2015

References