## Abstract

The effect of spatial aperture variations on the thermal performance of discretely fractured geothermal reservoirs was investigated using finite element method solutions to the convective heat transfer in the fracture coupled with a boundary integral description of conductive heat transfer in the rock. The dipolar flow was generated by a source 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 realizations was evaluated. Fractures with self-affine aperture fields where long range correlations 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.

## Background

Variations in fluid velocity on various flow paths in discretely fractured geothermal reservoirs like those found in Enhanced/Engineered Geothermal Systems (EGS) have a significant effect on the thermal performance of a reservoir. One source of variations in fluid flow is from spatial aperture variations in the fractures. EGS reservoirs are developed 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 pressurization 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 thermal 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 bulk rock. Because of the low thermal diffusivity of the host rock (of order \(10^{-6}\) m^{2}/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-affine 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 reservoir’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 streamlines where the fluid velocity is faster near the injection and production wells. For simplicity, effects like temperature dependent properties and buoyancy were not considered. We are interested in modeling reservoirs with different degrees of fracture surface variations, roughness/Hurst coefficients, wellbore spacings, and multiple fractures. Our focus is to determine how fracture aperture variations affect reservoir parameters like the production 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 provides a method to systematically generate variations in flow fluid rates and their pathway through the fracture system.

While little work has been done on the coupled convective-conductive heat transport 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 channeling 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) determined 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 distribution 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 (Auradou 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 surface 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 reservoir’s thermal performance. Jupe et al. (1995) attempted to model spatial aperture variations 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 variogram. 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 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^2/\alpha\), where *w* is the fracture aperture and \(\alpha\) is the thermal diffusivity of the rock. For an aperture value from 1 to 5 mm and a thermal diffusivity of \(1\times 10^{-6}\) m^{2}/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 compared 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 knowing 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 transport 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 temperatures will offer a better and more practical understanding of the effects fracture surface 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 \(\left\langle w \right\rangle\) to give

The spatial mean aperture \(\left\langle w \right\rangle\) is given by

where *A* is the fracture surface area and the variance or mean square of the aperture field is given by

The spatial correlation of an aperture field can be described by the variogram function \(\gamma\), which is twice the difference of the aperture variance and the two-position autocorrelation function:

where \(\mathbf {r'}\) is the displacement from point \(\mathbf {r}\) and \(\mathbf {r}\) is a point on the surface of the fracture. The autocorrelation function is a measure of the correlation of the aperture at point \(\mathbf {r}\) with that at point \(\mathbf {r}+\mathbf {r'}\). If the variogram function follows

it has a scaling behavior of \(\gamma (\lambda \mathbf {r'})=\lambda ^{2\zeta }\gamma (\mathbf {r'})\), where \(\lambda\) is a scaling parameter and the field is said to have a self-affine geometry governed by the Hurst or roughness coefficient, \(\zeta\) (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 materials 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 \(\zeta\) 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 wavelengths) are more important than small scale (large frequencies, short wavelengths) variations. The large scale variations become more strongly dominant when \(\zeta\) is increased. However, the large scale variations are the primary determinant of the aperture variations for any positive \(\zeta\). A self-affine aperture field is analogous to fractional Brownian motion, if one views the spatial displacement \(|\mathbf {r'}|\) as a time-like variable and the change of aperture with displacement as equivalent to the spatial step of the Brownian motion. With this analogy, \(\zeta <1/2\) corresponds to a subdiffusive or antipersistent process, \(\zeta = 1/2\) corresponds to a normal diffusive process, and \(1/2<\zeta <1\) corresponds to a super diffusive process with \(\zeta = 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 Transforms 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 2*R* 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 \(\sigma _h\) of the fracture, the aperture value *h* is rescaled.

Real reservoirs are not static; there will be dynamic changes to the aperture and permeability fields due to thermal, hydraulic, mechanical, and chemical (THMC) effects. In earlier works, these effects have been incorporated to model a reservoir’s thermal performance 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 viable 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 injection 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, precipitation, 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 sample 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

where *P* is the pressure, \(\mathbf {v}\) is the velocity, \(\rho _w\) is the density of water, and \(\mu\) is the dynamic viscosity. The conservation of mass can be expressed as

The summation on the right hand side refers to sources/sinks of mass, *w* is the local aperture or gap of the channel, and \(\left\langle \mathbf {v} \right\rangle\) is the fluid velocity averaged across the gap. The 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

and this velocity is a function of the position \(\mathbf {r}\) in the two-dimensional plane of the fracture. Combining Eqs. (7) and (8) yields

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 approximate lubrication analysis. As such, Eq. (9) is valid when \(\mathrm {Re}(w/L_{x,y}) \ll 1\) and \(w/L_{x,y} \ll 1\) (Deen 1998). \(L_{x,y}\) refers to the characteristic length scale of aperture variations in the *x* and *y* direction and may be considered a “wavelength” of aperture variations. Here, the Reynolds number is defined as \(\text{Re}\equiv (\rho _w U_c w)/\mu\), where \(U_c\) is the characteristic velocity. For a mass flow of 40 kg/s, a wellbore spacing *L* of 400 m, a width of 1 mm, and a characteristic velocity of \(U_c \sim \dot{m}/(\rho _w w L)\), then \(\text{Re}\sim 100\), well below the critical Reynolds number for turbulence transition. The lubrication criteria \(\text{Re} (w/L_{x,y}) \ll\) and \(w/L_{x,y} \ll 1\) will also be satisfied when \(\text{Re}\sim 100\) because \(L_{x,y} > 1\) m because large scale aperture variations dominate when the aperture fluctuations are self-affine. Eq. (9) has been extensively 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 pressure. 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 limited 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 element 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 reservoir. 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 for the flow, mass, and heat transport in a fracture, taking advantage of not having to discretize a three dimensional domain.

The fluid energy balance inside the fracture is

The term \(q_s\) refers to the surface heat flux and \(c_w\) and \(k_w\) are the heat capacity and thermal conductivity of water, respectively. By neglecting the transient heat storage in the liquid contained by the fracture and thermal conduction in the liquid (\(\mathrm {Pe}\equiv U_c L/\alpha _w \gg 1\)), Eq. (10) becomes

Unfortunately, the preceding equation cannot be solved independently because the surface heat flux is unknown and coupled to transient conduction of the rock mass surrounding the fracture. The surface heat flux is also related to the temperature of the system 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:

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 subscript *r* refers to the thermophysical properties of the rock while \(T_r\) is the original reservoir rock temperature. The Laplace transform of the fracture heat balance equation is simply

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. Equations (13) and (14) are solved simultaneously to obtain the surface heat flux and the temperature 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 \(\bar{f}(s)\) at different values of *s*.

The matrix \(\bar{G}_{ij}(x-x',y-y',z)\) is computationally intensive to construct because every entry must be evaluated. The computational cost is \(O(N^2)\) where *N* is the number of nodes. Computations can be greatly reduced by noticing that many of the entries in \(\bar{G}_{ij}\) are virtually zero and can be skipped in populating \(\bar{G}_{ij}\). For determining \(\bar{G}_{ij}\) for a given *s*, entries that correspond to evaluating points that are farther than \(20\sqrt{\alpha /s}\) are skipped where \(\alpha\) 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 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 \(\bar{G}_{ij}\). Additionally, the 1D simplification corresponds to a 400-fold decrease in simulation time compared to the fully populated \(\bar{G}_{ij}\) matrix. Fox et al. (2013) showed that 1D conduction is a reasonable approximation when \(\beta ^2 \ll 1\). Here, \(\beta\) is a dimensionless parameter describing the ratio of the extent of thermal depletion normal to the plane of the fracture to the extent of depletion in the plane of the fracture. For dipole flow, \(\beta\) can be defined as \(\beta = 2k_r L/(\dot{m} c_w)\), where *L* is the spacing between the injection and production well. For typical flow conditions of 40 kg/s and a wellbore separation of 400 m, \(\beta ^2 = 2 \times 10^{-4}\). Nonetheless, for completeness, we used the full 3D heat conduction in this study.

The dimensionless temperature in the reservoir is defined as \(\Theta \equiv (T(x,y)-T_w)/(T_r-T_w)\), where *T* is the temperature of the fracture surface at point (*x*, *y*), \(T_w\) is the temperature of the injected water and \(T_r\) is the initial rock temperature which remains constant in the far-field. The value of \(\Theta\) 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 \(\Theta =1\), there is no drawdown since the temperature in the fracture is still at the initial reservoir rock temperature. The variable \(\Theta _p\) refers to the dimensionless production temperature. If \(\beta ^2 \ll 1\) (1D conduction), then the dimensionless temperature is a function of the dimensionless wellbore spacing and dimensionless time:

where \(\tau \equiv L^4 k_r \rho _r c_r/(\dot{m} c_w)^2\). The dimensionless wellbore spacing is constrained to range from zero to one. The results for one mass flow rate and a set of times for a given dimensionless wellbore spacing are the same as for a different mass flow rate and a set of extraction 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, \(\sigma _h/\left\langle w \right\rangle\) and \(\zeta\). For this study, we have used a radius *R* of 500 m and an injected mass flow \(\dot{m}\) of 40 kg/s. The thermophysical property values used were \(k_r=2.4\) W/(m °C), \(\rho _r=2300\) kg/m^{3}, \(c_r=1000\) J/(kg °C), and \(c_w=4184\) J/(kg °C). The value of \(\tau\) is 2.1 \(\times\;10^7\) s or 0.16 years for a wellbore spacing of 400 m and for the previously 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 modeling provides the statistical distribution of the dimensionless production temperatures as a function of wellbore spacing, duration of extraction, the variance of fracture aperture 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 \(\Theta _{p,0}\) refers 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 \(\Theta _{p,0}=0.8\) is a realistic value for the drawdown criteria in many practical situations.

For \(\sigma _h/\left\langle w \right\rangle\), \(\zeta\), and *L*/(2*R*), we chose a range of possible values for these parameters to model their impacts on thermal performance. The range of \(\sigma _h/\left\langle w \right\rangle\) was chosen to be from 0.015 to 0.35. Granite core samples have been shown to have a value of \(\sigma _h/\left\langle w \right\rangle\) of 0.42 (Hakami and Larsson 1996) but also low values around \(5\times 10^{-4}\) (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 \(\sigma _h/\left\langle w \right\rangle = 0.35\), it becomes difficult to generate fractures that do not have zero or negative aperture values. Thus, the aperture realizations 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 coefficient 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*/(2*R*) 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 \(\sigma _h/\left\langle w \right\rangle\) with *L*/(2*R*) and \(\zeta\) kept constant at 0.4 and 0.8, respectively. In a second set of simulations, we repeated the process but now varying the Hurst coefficient and keeping \(L/(2R)=0.4\) and \(\sigma _h/\left\langle w \right\rangle =0.25\). Finally, in a third set, the wellbore spacing was varied and \(\sigma _h/\left\langle w \right\rangle\) and \(\zeta\) are set to 0.25, and 0.8, respectively. The statistical values of the thermal performance of most interest are the mean dimensionless production temperature \(\bar{\Theta }_p\) of the ensemble, the standard deviation of production values \(\sigma _{\Theta }\), 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 generate 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

For a range of \(\sigma _h/\left\langle w \right\rangle\) from 0.015 to 0.35, five thousand fracture realizations were generated for each value of \(\sigma _h/\left\langle w \right\rangle\). For each five thousand fracture realization data set, the values of \(\bar{\Theta }_p\), \(\sigma _{\Theta }\), and *F* are reported based on the thermal performance for all realizations at the time when the dimensionless production temperature reached 0.8 in the 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 dimensionless heat transfer area accessed at a certain time is defined as \(\int (1-\Theta )\ \mathrm {d}A/A=1-\left\langle \Theta \right\rangle\). Thus, the mean fracture surface temperature \(\left\langle \Theta \right\rangle\) is a measure of the thermally accessed area of the fracture at a given extraction time. Comparing \(\left\langle \Theta \right\rangle\) for different fracture aperture distributions, a fracture with a lower value of \(\left\langle \Theta \right\rangle\) 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 \(\left\langle \Theta \right\rangle\) is 0.816 when \(\Theta _{p,0} = 0.8\). For a fracture with smaller apertures in the region between the wells, the fluid is pushed to higher aperture regions in the periphery and the effective heat transport area increases, \(\left\langle \Theta \right\rangle =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 determining whether the reservoir will exhibit poor or good thermal performance.

Figure 2 shows the distribution of thermal performance when \(\sigma _h/\left\langle w \right\rangle =0.25\). Although a Gaussian curve does not fit the histogram perfectly, it does an adequate job of describing the behavior. A Gaussian distribution of \(\Theta _p\) is expected for a Gaussian aperture fluctuations with small amplitude variations. Inserting \(\left\langle w \right\rangle + h\) into the governing equations for fluid velocity and fracture surface temperature and considering a regular perturbation expansion for small \(\epsilon = \sigma _h/\left\langle w \right\rangle\), one obtains \(\Theta _p=\Theta _0+ \epsilon \Theta _1 + \epsilon ^2 \Theta _2 + O(\epsilon ^3)\) for *w*, where \(\Theta _0\) is the production temperature for the constant aperture case and \(\Theta _1\) is a linear functional of the field \(h(\mathbf {r})/\sigma _h\). A linear functional of a Gaussian field is a Gaussian variable. 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 aperture variations affects the mean and standard deviation of the production temperature. As \(\sigma _h/\left\langle w \right\rangle\) 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 \(\sigma _h/\left\langle w \right\rangle\).

The linear behavior of Fig. 4 is expected due to the before mentioned linear behavior of \(\Theta _p\) to small variations of *h*. In particular, \(\Theta _0\) has zero standard deviation, while \(\Theta _1\) has an *O*(1) standard deviation so that the standard deviation of \(\Theta _p \approx \Theta _0+ \epsilon \Theta _1\) is proportional to \(\sigma _h\). The quadratic relationship between the mean production temperature and the aperture standard deviation is also expected owing to the fact that \(\bar{\Theta }_1\) must be zero based on linearity. Thus, \(\bar{\Theta }-\Theta _0 = (\sigma _h/\left\langle w \right\rangle )^2 \bar{\Theta }_2\), where \(\bar{\Theta }_2\) is a constant of proportionality.

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 \(\Theta _p\) reaches 0.8. When \(\sigma _h/\left\langle w \right\rangle =0.1\), which resulted in \(\sigma _{\Theta }=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 \(\sigma _h/\left\langle w \right\rangle =0.35\), corresponding to \(\sigma _{\Theta }=0.092\), the results 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 \(\sigma _h/\left\langle w \right\rangle\); when \(\sigma _h/\left\langle w \right\rangle =0.05\), the 25 and 75 % quantiles of reservoir life are closer together and equal to 0.91 and 1.08, respectively.

For a normally distributed production temperature, the probability density function is equal to

From Fig. 2, it is reasonable to model the production temperature using a normal (Gaussian) distribution. The fraction of fractures that exhibit enhanced thermal performance is equal to one minus the integral of Eq. (16) with respect to \(\Theta _p\) form 0 to \(\Theta _{p,0}\). The result is

As discussed earlier, the mean and standard deviation of the production temperature distribution have a quadratic and a linear behavior, respectively:

and

Using Eqs. (18) and (19), the fraction of fractures that exhibit enhanced thermal performance compared to the base case can be expressed in terms of the standard deviation of the aperture fluctuations,

where *K* is the ratio of \(C_1/C_2\). Since the error function is linear for small values of \(K \sigma _h/(\sqrt{2} \left\langle w \right\rangle )\), Eq. (20) can be approximated as

Figure 5 displays the linear decrease of *F* as a function of \(\sigma _h/\left\langle w \right\rangle\). 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 \(\Theta _p\) is Gaussian. The theoretical and simulated curves only deviate about 0.015 at the end of the range investigated. The fraction 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 different 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 production temperature can be modeled as a Gaussian distribution. The fraction of fracture 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 \(\sigma _h/\left\langle w \right\rangle\) 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 \(\zeta\) on thermal performance, simulations were performed for ten sets of five thousand fractures. Each set had a different value of \(\zeta\) in the range from 0.1 to 1 and \(\sigma _h/\left\langle w \right\rangle\) was kept constant at 0.25. Figure 6 displays 1D fracture aperture profiles for various values of \(\zeta\). Increasing \(\zeta\) yields aperture fields with strong large scale correlations. In other words, the profile for \(\zeta =1\) is smoother than the profile for \(\zeta =0.1\), which is characterized by many “peaks” and “valleys”. Table 1 summarizes the statistical results of the ensembles for each value of \(\zeta\), with a 95 % confidence interval derived by resampling based on the bootstrapping method (Efron and Tibshirani 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 coefficient 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 variance 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 \(\sigma _{\Theta }\) changes appreciably, ranging from 0.063 to 0.073 with increasing \(\zeta\). If the aperture field is more correlated (larger value of \(\zeta\)), the 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 \(\zeta\). The standard deviation of the production temperature does not change much for \(\zeta >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*/(2*R*). As *L*/(2*R*) 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*/(2*R*) of 0.2, 0.4, 0.5, 0.6, and 0.8, the thermal performance of the five thousand fracture simulations for each value of *L*/(2*R*) were all performed for an end time when the respective *L*/(2*R*) base case production temperature reaches 0.8. To reach the \(\Theta _p=0.8\), the dimensionless time \(t/\tau\) for *L*/(2*R*) 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/\tau\) to reach the drawdown criteria drops for increasing values of wellbore spacing because \(\tau\) is proportional to \(L^4\). The values of \(\zeta\) and \(\sigma _h/\left\langle w \right\rangle\) were kept constant at 0.8 and 0.25, respectively.

Boxplots for the distribution of thermal performance for the different values of *L*/(2*R*) are plotted in Fig. 7. The boxplots show that the distribution of dimensionless production 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 temperature is greatest for this value. The mean value of the distribution and *F* decreases as *L*/(2*R*) 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 \(\sigma _{\Theta }\). 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 temperature \(\left\langle \Theta \right\rangle\) is a measure of the accessed area for heat transfer at a point in time; smaller values of \(\left\langle \Theta \right\rangle\) mean more area has been accessed. The base case values for \(\left\langle \Theta \right\rangle\) 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.

### 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 \ll 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 \(\dot{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 baseflow 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 distribution of flow in the reservoir.

Since the pressure drop \(\Delta P\) is linear with the flow rate, \(\Delta P \propto \dot{m}\), the flow split fractions for fracture A and B are

where \(I_i \equiv \Delta P/\dot{m}_i\) and represents the flow impedance of fracture *i*. The system can be 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 modify 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 \(\phi\) as \(\mathrm {max}(\phi _A,\phi _B)\), which ranges from 0.5 to 1.

The hydraulic aperture is another way to view the varying pressure drop for each fracture. In terms of the pressure drop of a fracture \(\Delta P\), the hydraulic aperture is defined as

where \(\Delta P_0\) is the reference pressure drop in a smooth fracture for an aperture value of \(w_0\) 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 studied by determining the thermal performance of an ensemble of fracture realizations. The wellbore spacing and the Hurst coefficient were kept constant and the standard deviation of the aperture was varied over the same range considered before. Finally, we examined 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 \(\zeta =0.8\) were generated and their thermal performance was simulated. Reassigning the pairing had an insignificant effect on the results presented. Figures 8 and 9 show the histogram of the production temperature and the flow split fraction when \(\sigma _h/\left\langle w \right\rangle =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 \(\phi\) values, corresponding to more unequal flow splits, are likely to result in lower production temperatures. Also, plotted in the figure is the production temperature as a function of the flow split if the fractures had no surface variations. In the absence of aperture variations, \(\Theta _p\) is always monotonically decreasing with respect to \(\phi\). 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 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 thermal performance has a linear relationship with the uncertainty in \(\phi\).

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 \(\sigma _{\Theta }\) and \(\bar{\Theta }_p\). Compared to the single fracture case, the values of both \(\sigma _{\Theta }\) and \(\bar{\Theta }_p\) are smaller. The curve is no longer linear since *K* is now equal to 3.00 and as a result, the linear approximation of the error function breaks down.

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(\phi )\). 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 \(\Theta _p > g(\phi )\), 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 \(\phi\). However, \(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 performance of the two fracture system is the aperture variation and flow channeling within 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 individual fracture. In general, a fracture with a small aperture region in between the injection 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_H\) and \(\Theta _p\) for the single fracture is \(-0.65\) when \(L/(2R)=0.4\), \(\zeta =0.8\), and \(\sigma _h/\left\langle w \right\rangle =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 compared with single fracture reservoirs. However, the negative impact on thermal performance 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, 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\), \(\sigma _h/\left\langle w \right\rangle =0.25\), and \(\zeta = 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 \(\sigma _h/\left\langle w \right\rangle =0.015-0.35\), *F* monotonically decreases from 0.49 to 0.074. The three fracture 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 fracture 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 fracture 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 distribution of the two fracture system when \(L/(2R)=0.4\), \(\zeta =0.8\), and \(\sigma _h/\left\langle w \right\rangle =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 flow split represents the value of \(\phi\) that optimizes the production temperature. In practice, 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 performance. 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 temperature 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 \(\Delta \Theta _p\) versus the flow split fraction. For the 50/50 flow split scenario, the results for \(\Delta \Theta _p\) versus \(\phi\) 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 \(\Theta _p\) 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 improvements of applying flow control increases; the covariance of \(\phi\) and \(\Delta \Theta _p\) 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 \(\sigma _h/\left\langle w \right\rangle =0.05\), 75 % of the cases only had a \(\Delta \Theta _p\) of less than 0.002 when the optimal value of \(\phi\) was used.

Compared to the *F* value for the one fracture reservoir, applying flow control for the range of \(\sigma _h/\left\langle w \right\rangle\) 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 \(\sigma _h/\left\langle w \right\rangle >0.25\) that *F* of the equal flow split is appreciably lower than the one fracture reservoir.

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 thermal 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 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 connected 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 that the dimensionless standard deviation of the aperture \(\sigma _h/\left\langle w \right\rangle\) is more important in controlling thermal performance than the Hurst coefficient \(\zeta\). Specifically, both the uncertainty in thermal performance, \(\sigma _{\Theta }\), 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 insignificant when \(\sigma _h/\left\langle w \right\rangle < 0.1\). For the investigated range of \(\sigma _h/\left\langle w \right\rangle\) from 0.015 to 0.35, 34 to 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 aperture 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 rectangular 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 discretely 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 deterioration in thermal performance arises because fractures with a larger flow impedance usually 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 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 aperture field stayed constant with time. Real reservoirs, however, will experience thermal, 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 perpendicular 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 currently 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 reservoir 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 (Castagna 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 aperture, 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.

## 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 fractures: 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 penetrating 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 fractured 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 stochastic 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 extraction 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.

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 influence 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 geological 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 numerical 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. Effect 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 thermoelastic 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.

## Authors’ contributions

This study was performed by DBF with advice on the planning, design, and geothermal application of the model provided 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.

### Acknowledgements

The authors would like to thank the National Science Foundation IGERT program (contract 0966045), the U.S. Department 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.

## Author information

### Affiliations

### Corresponding author

## Rights and permissions

**Open Access** 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.

## About this article

### Cite this article

Fox, D.B., Koch, D.L. & Tester, J.W. The effect of spatial aperture variations on the thermal performance of discretely fractured geothermal reservoirs.
*Geotherm Energy* **3, **21 (2015). https://doi.org/10.1186/s40517-015-0039-z

Received:

Accepted:

Published:

### 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