Numerical Investigations on the Impact of Fracture Aperture Anisotropy on EGS Thermal Performance

The overall goal of this study was to investigate how fracture aperture anisotropy can impact flow and heat transport , and to demonstrate ways Enhanced Geothermal Systems can be harnessed to optimize thermal performance. To achieve the goal of this study, a systematic fracture characterization approach was used, and numerical simulation models were used to study the physical processes that govern the interaction between the fluid and the rock during heat extraction from Enhanced Geothermal Systems. It was demonstrated in this study that the flow-wetted surface area had a direct and significant contribution to the amount of heat extracted. For the lab-scale fractures, the Joint Roughness Coefficient (JRC) confirmed geometric anisotropy of the fracture aperture and was seen to have a direct correlation with the flow contact area. The lower the difference in JRC values between the perp endicular and parallel flow configurations, the more flow contact area is expected in the perpendicular flow direction, which will lead to more heat extracted from the rock. From the variogram model parameters, it was deduced that high geometric anisotropy results in high differences in thermal drawdown. The thermal performance appeared to be better in the perpendicular flow configuration with a ratio of 70:30 for the lab-scale fractures. For the field-scale fractures, it was seen that most of the synthetically generated fracture aperture distributions with a geometric anisotropy ratio of 2 had Hurst exponents of fracture surface aperture variability found in nature. For all the fracture aper ture distributions analyzed for the field scale, the perpendicular flow configurat ion resulted in better thermal performance than the parallel flow configuration with a ratio of 64:36. Furthermore, for the geometric anisotropy ratio of 2, the ratio was 70:30. The results of this study suggest that placing an injector well in the direct ion perpendicular to shear or slip of an enhanced geothermal system may result in favorable thermal performance over a parallel flow configuration.


INTRODUCTION
Natural fractures are typical features in many rocks, occurring at various length scales. It has been documented that these fractures have surfaces with spatial variations in aperture that lead to tortuous flow pathways (Brown, 1987;Tsang & Neretnieks, 1998). Understanding the spatial variations of fracture apertures and their impact on fluid flow is relevant for systems that rely on physical processes dependent on flow, such as unconventional reservoirs and geothermal reservoirs. Spatial aperture variation is often referred to as "roughness" in the literature. Brown (1989) investigated the relationship between fractures' hydraulic properties and electrical properties with spatial variations. The results indicated that fluid flow and electric current movement through fractures with spatial variations in aperture are primarily dependent on the area of contact between the fluid and the rock surface, while the details of the fracture aperture variability were of secondary importance. The fracture aperture variability analyzed considered the fracture aperture distribution to be isotropic.
Anisotropic roughness is often observed on natural fracture surfaces, and surface roughness anisotropy may also manifest in the fracture aperture (Thompson and Brown, 1991). Candela et al. (2012) examined the topographic aperture variability measurements of five exhumed faults, and 13 surface earthquake ruptures with scales ranging from 50 μm to 50 km. They observed that fault aperture variability is scale-dependent, with an anisotropic self-affine behavior. They determined the Hurst exponent, H, which describes the roughness of the fracture surface, to be 0.58 ± 0.07 in the slip direction and 0.81 ± 0.04 in the direction perpendicular to slip (Schmittbuhl et al., 2008, Candela et al., 2012. The impact of fracture aperture anisotropy on flow properties has been studied. Thompson and Brown (1991) investigated the role of anisotropic surface aperture variability on fluid flow, solute transport, and electrical current flow in fractures. They determined that aperture variability oriented parallel to the primary flow direction enhanced fluid and solute transport rates, while aperture variability oriented transverse to the flow direction inhibited flow rates and delayed solute movement through the fracture. Co (2017) analyzed flow properties of lab-scale fractures. The study by Co (2017) determined that sheared fractures exhibited geometric and permeability anisotropy. Using Sequential Gaussian Simulation and the variogram from the lab-scale fractures, Co (2017) generated several artificial aperture distributions to further investigate the lab-scale fractures' geometric and permeability anisotropy. Flow configurations were assigned based on the lateral direction of the shear offset for the sheared fractures. Thus, the perpendicular flow configuration had a pressure drop perpendicular to the lateral shear offset direction. In contrast, the parallel flow configuration had a pressure drop parallel to the lateral shear offset direction. From the study by Co (2017), it was determined that for the perpendicular flow configuration, 97% of flow occurred in 26% of the fracture area, while for the parallel flow configuration, 97% of the flow occurred in 15% of the fracture area, implying that there was lower contact area in the parallel flow configuration.
Enhanced Geothermal Systems rely on heat extraction by the physical process of fluid flow. Colder fluid is injected into hot rock with flow predominantly through fractures. Heat from the rock surrounding the fractures is transferred to the production well(s) through the flowing fluid that encounters the rock; hence only the portion of the fracture in contact with the flowing fluid provides effective heat exchange surface area. The more heat exchange surface area contacted by the flowing fluid, the more effective is the heat ext raction.
Several studies have been carried out to understand the coupled flow and heat process in Enhanced Geothermal Systems. On the laboratory scale, such studies include the works by Zhao & Tso (1993), Huang et al. (2016), Bai et al. (2017Bai et al. ( ), M a et al. (2019, and Zhang et al. (2020). On the commercial scale, such studies include the works by Neuville et al. (2010), Hawkins &Becker (2012), andFox et al. (2015). These studies demonstrated that the geometry of the fracture surface has significant implications for heat transfer through the fracture. However, these studies did not consider the heat transfer implications of fracture aperture variability anisotropy. Gao et al. (2021) studied the influence of flow direction on the heat transfer characteristics applicable to a single fracture in granite. They established a single fracture heat transfer model with a random geometry profile. Four cases with fracture profiles and varying angles between flow directions were set up to simulate and explore the heat transfer performance of distilled water through fractures. They used the parameter α of values 0 o , 30 o , 60 o , and 90 o for the surface flow direction. 0 o was parallel to the direction of flow, while 90 o was perpendicular to flow direction. From their study, they deduced that the fracture surface with α = 0 o resulted in more cooling than the fracture with α = 90 o .
The fracture surface used in the study by Gao et al. (2021) was a pseudo-3D fracture surface that does not account for the anisotropy of the fracture aperture variability. The fracture aperture distribution thus may not represent the actual nature of real fractures. Also, only four fracture surface distributions were investigated. Hence, in this study, the impact of anisotropy of fracture aperture variability on heat transfer was investigated using fracture surface aperture variability derived from real laboratory fractures and included an additional 210 artificially generated fracture aperture distributions representative of real fractures.

Model Description
The system modeled is a hypothetical EGS doublet consisting of a single injector/producer well-pair circulating the working fluid through a single fracture contained within hot, impermeable, granitic rock. Two model scales were studied. On the laboratory scale, the fracture was 50 mm x 50 mm in the horizontal plane embedded within a relatively impermeable bulk rock matrix of height 100 mm. Hypothetical horizontal wells, one injection and one production, were placed at the edges of the fracture. The numerical model consists of a 50 by 50 by 45 grid. In the horizontal x and y directions, the individual cells were of uniform length of 1 mm, while in the vertical z direction, the thicknesses were very fine around the fracture and become coarse away from the fract ure. Figure 1 shows a snapshot of the reservoir simulation domain. On the field scale, the fracture was horizontal, measuring 1000 m x 1000 m at a depth of 1295.5 m below ground surface, and was embedded within the relatively impermeable bulk rock matrix. Horizontal wells, one injection, and one production, were placed at the edges of the fracture. The numerical model is a 50 by 50 by 45 grid. In the X and Y directions, the individual cells are of uniform length of 20 m, while in the Z direction, the thicknesses were very fine around the fracture and become coarser away from the fracture. Figure 2 shows a snapshot of the field-scale reservoir simulation domain.  The coupled flow and heat transport mechanism was modeled with a three-dimensional compositional numerical simulator -ECLIPSE. ECLIPSE is a finite-difference simulator and was run in the fully implicit mode, using Cartesian block-center geometry in three dimensions for flow and heat transport. The simulator has been verified for geothermal applications by Stacey and Williams (2017) and Okoroafor and Horne (2018).
The initial reservoir conditions for the lab-scale simulation were 200 °C and 0.1 M Pa, while for the field-scale simulations, the initial conditions were 200 °C and 10 M Pa. Table 1 shows the rock and fluid properties and relevant parameters used in the model. Viscosity and density were treated as temperature-dependent. Okoroafor and Horne

Model Assumptions
The following assumptions were considered in setting up the system to be modeled with water as the working fluid: -


The fluid circulating throughout the system is single-phase and remains in the liquid state throughout the simulation.
 Fluid flow is in the laminar regime with Reynold's number low enough to apply Darcy's law.


The impact of chemical dissolution/deposition, thermal stresses and changes in aperture were ignored.

Fracture Characterization
The fracture is treated as a porous medium with porosity set as 0.99 while the heterogenous permeability is defined by the local cubic law for a fracture with spatial variations (Oron & Berkowitz , 1998), which is represented by Equation 1. (1) where kf, i, j, b, and H are the effective permeability, grid number in the x-direction, grid number in the y -direction, local fracture aperture, and thickness of the fracture grid element, respectively.

Lab Scale Aperture Distribution Determination
The work by Ishibashi et al. (2012) provides nonuniform aperture fields for sheared fractures. Co (2017) derived variogram model parameters of the heterogeneous aperture field from Ishibashi et al. (2012), which were subsequently used to generate artificial aperture fields. This was done using Sequential Gaussian simulation (SGSIM ) with the Stanford Geostatistical M odeling Software (SGeM S) (Remy et al., 2009). A full discussion on the SGSIM method and variogram modeling can be found in Goovaerts (1997). Figure 3 shows the selected aperture distribution from Ishibashi et al. (2012) used to generate the artificial aperture distributions. The figure also shows directions to indicate a parallel flow configuration and a perpendicular flow configuration.

Field Scale Aperture Distribution Determination
The 75 mm x 50 mm lab-scale aperture distribution from Ishibashi et al. (2012) was the starting point to generate the field-scale aperture distribution. First, the maximum aperture for the field scale was determined using the scaling correlation by Olsen (2003), which Co (2017) fitted to give Equation 2. (2) Using the calculated maximum aperture, the aperture distribution was adjusted to match the maximum aperture for the field scale. A variogram model was also determined from the lab-scale data. Subsequently, synthetic fracture aperture maps were generated using Sequential Gaussian Simulation (SGSIM ). The main input data sets for the SGSIM runs were the adjusted aperture distribution and the aperture variogram model. Then the aperture distributions were checked to determine if their Hurst exponents were in the range found generally in nature. In this study, Hurst exponents used for selecting the aperture distribution were 0.6 ± 0.1 in the slip direction and 0.85 ± 0.1 in the direction perpendicular to slip. Once an aperture distribution met these criteria, its permeability was estimated for input into the thermohydraulic model.

Joint Roughness Coefficient (JRC), Hurst Exponent and Variogram Modeling
The joint roughness coefficient (JRC) was estimated using a fracture surface analysis toolbox in M ATLAB by Heinze et al. (2021), which automatically compares the estimated JRC against the 10 standard roughness profiles defined in the work by Barton and Choubey (1977). The higher the JRC, the higher the spatial aperture variability of the surface will appear. Heinze et al. (2021) codes were used to estimate the Hurst exponent and fractal dimensions of two-dimensional and three-dimensional surfaces with aperture variability. The method calculates the Hurst exponent and the fractal dimension based on the work by Schmittbuhl et al. (2008).
The variogram is a measure of spatial dissimilarity. Variogram modeling involves fitting a variogram model to a set of data using models known before modeling that meet a positive-definiteness criterion to ensure the interpolation method done on the dataset results in a unique solution. The following are examples of positive-definite theoretical semivariogram models commonly used: the nugget effect, spherical, exponential, Gaussian, and power models (Goovaerts, 1997). In this study, the exponential semivariogram model was used to fit the analog lab-scale fracture aperture data to provide an experimental variogram that could be used for further analysis. To account for variable spatial continuity in different directions, which is the geometric anisotropy, separate empirical and model variograms can be estimated for different directions in the data set.

Impact of fracture anisotropy on lab-scale fractures
The thermohydraulic model was run with the lab-scale fracture characterized by 100 artificially generated fracture aperture distributions. Figure 4 is a plot of the difference in thermal drawdown between the perpendicular flow direction and the parallel flow direction for all 100 simulations. All plots above the zero horizontal line indicate that the temperature measured at the extraction end of the fracture was higher in the perpendicular direction than in the parallel direction. Similarly, all plots below the zero line indicate that the temperature measured at the extraction end of the fracture was higher for the parallel flow configuration resulting in lower thermal drawdown in that flow configuration. From Figure 4, it can be seen that the highest temperature differences occurred between 0.6 and 0.8 minutes into the simulation. Hence a histogram of the temperature differences at 0.72 minutes was plotted in Figure 5. The data in Figure 5 indicates that most of the temperature differences are above 0, and the temperature difference values with the highest frequency lie between 1.25 °C and 2.5 °C. The data were sorted, and a count was performed to quantify the percentage of the artificially generated fracture aperture distributions that resulted in the perpendicular flow configuration with lower thermal drawdown than the parallel flow configuration. Figure 6a shows the sorted data at 0.72 minutes of simulation, while Figure 6b shows the sorted data at 10 minutes. At 0.72 minutes of simulation, the percentage of fracture aperture distributions that resulted in the perpendicular flow configuration having lower thermal drawdown (and higher temperature at the extraction end of the fracture) was 70 %. By 10 minutes of simulation, the value was 71 %.
These results indicate that, on average, there is a 70 % chance that the perpendicular flow configuration will result in a lower thermal drawdown than the parallel flow direction.  Fig. 6a) and 10 minutes (Fig. 6b) into the simulation.
The temperature difference between the perpendicular and parallel flow directions for the 100 artificially generated fracture aperture distributions at 0.72 minutes was plotted against the difference in the fracture area contacted by the flowing fluid. The plot is shown in Figure 7. The data in the p lot indicates a linear relationship between the temperature difference between the flow directions and the difference in the fracture surface area contacted by fluid. Thus, the higher the difference in the flow wetted surface area, the higher the difference in temperature between the flow configurations. A few aperture distributions were selected for further analysis to understand more the impact of anisotropy on the thermal performance of the lab-scale fractures. Figure 8 shows the temperature difference plot of the chosen aperture distributions, while Table 2 lists the selected aperture distributions with remarks to support why they were chosen. These aperture distributions were analyzed based on their flow and temperature maps at the fracture, contact area, joint roughness coefficient (JRC), and anisotropy of the ranges from variogram modeling. The results of all these parameters are summarized in Table 3. The objective was to understand if specific characteristics of the fracture surface favored heat transport for a given flow configuration.  Figures 9 and 10 show the flow and temperature maps at the fracture surface for Aperture 11 and Aperture 63 at 10 minutes into the simulation. Table 3 summarizes the percentage of the fracture area covered during flowing conditions for all the selected aperture distributions.
Aperture 11 had the highest temperature difference recorded from the simulations. From the flow map of Figure 9, flow is seen to be more channelized in the parallel flow configuration (Figure 9a) than in the perpendicular flow configuration (Figure 9b). As a result, less area (5.84 %) is contacted in the parallel flow configuration than in the perpendicular flow configuration (Table 3). From the temperature map, the less cooled area is more prominent in the parallel flow configuration (Figure 9c) compared to the perpendicular flow configuration, contributing to why more heat is extracted in the perpendicular flow configuration.   Aperture 63 had the highest temperature difference in favor of the parallel flow configuration to the perpendicular flow configuration with the parallel flow configuration having a higher temperature of about 12.39 °C by 0.72 minutes into the simulation. The flow map and temperature map at the fracture is shown in Figure 10. A qualitative inspection of the flow map shows more flow areas in the parallel flow configuration than in the perpendicular flow configuration (Figures 10a and 10b). From Figure 10c, more of the fracture surface area is cooled in the parallel flow configuration than the perpendicular flow configuration (Figure 10d). There is evidence of an area in the latter temperature map that has not been adequately cooled relative to the parallel flow configuration. In analyzing the JRCs alongside the p ercentage of flow contact area, it was observed that the JRCs in the parallel flow direction were higher than the JRCs in the perpendicular flow configuration. This information gives an indication of anisotropy of the fractured surface and agrees with studies that suggest that there is anisotropy in the aperture variability of fractured surfaces (Thompson and Brown, 1991;Candela et al., 2012;and Co, 2017).
For the aperture distributions analyzed, a low difference in the JRC value between the perpendicular flow configuration and the parallel flow configuration correlates with a high difference in the percentage of flow contact area. Thus, the lower the difference in JRC values, the more flow contact area expected in the perpendicular flow direction, which will in turn lead to more heat extracted from the rock.
Geometric anisotropy is said to exist when the ranges of the variogram model vary as a function of direction. From the variogram modeling, there were differences between the sills and the ranges across the different aperture distributions (Table 3). It was also observed that Aperture 11 and Aperture 63, which had the highest differences in temperatures betw een the flow configurations, also had the highest geometric anisotropies (1.54 and 1.63 respectively) and had lower sills compared to the other aperture distributions. On the other hand, Aperture 35 and Aperture 56, which had low temperature differences for the flow configurations (including changing from favoring the perpendicular flow configuration to favoring the parallel flow configuration), had low anisotropy. Hence it can be deduced that high geometric anisotropy results in high differences in thermal drawdown. Table 4 shows the categories of aperture distribution generated for the field scale with their corresponding geometric anisotropy ratios and ranges. Figure 11 is a boxplot showing the ranges of the Hurst exponents for each of the geometric anisotropy ratios for the perpendicular and parallel directions. The cyan bar shows the range of = 0.6 measured from real faults and fractures while the magenta bar shows the range of = 0.8 measured from real faults and fractures. From Figure 11 it can be seen that the aperture distributions with geometric anisotropy ratio of 2 fall in the range of Hurst exponents for real faults and fractures, and hence the results will be presented here. Figure 12a shows the sorted temperature difference between the perpendicular and parallel flow configurations one year into the simulation. Figure 12b shows the sorted data five years into the simulation. At one year of simulation, the percent age of fracture aperture distributions that resulted in the perpendicular flow configuration having lower thermal drawdown was about 64.5 %. The value was approximately 56.3 % five years into the simulation. These results indicate that there might be differences in which flow configuration is better at different times in the life of an Enhanced Geothermal System.   Fig. 12a) and five years (Fig. 12b) into the simulation. The difference is temperature in the perpendicular direction minus temperature in the parallel direction measured at the extraction well.

Impact of fracture anisotropy on field-scale fractures
For the geometric anisotropy ratio of 2, the y -direction range values of 50 m, 60 m, 75 m, and 100 m correspond to the x-direction range of 100 m, 120 m, 150 m, and 200 m. Figure 13 shows the sorted count of the temperature difference between the perpendicular and parallel flow directions for the different cases with geometric anisotropic ratio of 2 taken at one year into the simulation. This geometric anisotropy ratio had a total of 47 aperture distributions. About 70.2% of the aperture distributions favored the perpendicular flow configuration for this geometric anisotropy ratio. Recall from Figure 12 that the geometric anisotropy ratio of actual fractures was between 2 to 3. Thus, real fractures may have a 70% chance of better thermal performance in the perpendicular flow configuration than the parallel flow configuration. Two aperture distributions were selected for further analysis. The first had a temperature difference of -7.8 °C between the perpendicular flow configuration and the parallel flow configuration one year into the simulation. Henceforth, it will be referred to as Aperture 1. The second had a temperature difference of 26.2 °C between the perpendicular flow configuration and the parallel flow configuration one year into the simulation. Henceforth, it will be referred to as Aperture 2. Both aperture distributions had the lowest and highest differences, respectively, between the temperatures in the perpendicular and parallel flow configurat ions. Figure 14 shows the flow maps and temperature maps for the parallel and perpendicular flow configurations of Aperture 1. From Figure 14a, the flow in the parallel direction seems to flow along several striations on the fracture. However, the perpendicular flow configuration (Figure 14b) has flow predominantly in one preferential path and fewer fluid flow paths than the parallel flow direction. The estimated flow wetted surface was 34.4% for the parallel flow configuration and 29.4% for the perpendicular flow direction. Figures 14c and 14d show the temperature maps for the two flow configurations of Aperture 1. There is evidence of thermal sweep in near-parallel paths across the fracture plane in the parallel flow configuration (Figure 14c). Several fluid flow paths in the parallel flow configuration allow for more heat extraction from the system than the perpendicular flow configuration, where the channelized flow is limited to one dominant flow path and few alternative fluid flow paths. Thus, the pres ence of interconnected flow paths leading to more flow wetted surface area contributed to making the parallel flow configuration yield higher temperatures measured at the producer than the perpendicular flow configuration.
The results for the aperture distribution with the highest difference in temperatures at one year into the simulation between perpendicular and flow configuration are shown in Figure 15 for Aperture 2. In both flow configurations, the flow was highly channelized within preferential paths (Figures 15a and 15b), leading to large areas of low heat extraction (Figures 15c and 15d).  It was of interest to know why Aperture 2 led to highly channelized flow. Figure 16 shows the aperture distribution for Aperture 2 across the fracture plane with indications for perpendicular flow configuration and parallel flow configuration relative to t he direction of slip/shear. The bottom part of the aperture distribution has areas of non-zero aperture that are not connected. The apertures at the top part of the aperture distribution were connected, hence the channelized flow and heat transfer behavior seen in Figure 15.

CONCLUS IONS
The effect of fracture geometric anisotropy on thermal performance has been investigated by examining heat transport in parallel and perpendicular flow configurations defined by applying a pressure gradient relative to the direction of the fracture shear offset. A thermohydraulic model was used to represent the lab and field scales' behavior.
Of the 100 artificially generated fracture aperture distributions used in the lab-scale study, 70 % had higher temperatures in the perpendicular flow configuration than the parallel flow configuration. It was also observed from this study that the flow wetted surface area had a direct and significant contribution to the amount of heat extracted, i.e., the larger the surface area in contact with the flowing fluid, the more heat is extracted.
Geometric anisotropy was observed in the datasets joint roughness coefficient (JRC) and variogram modeling parameters. The JRC was seen to directly correlate with the flow contact area. The difference in JRC between the perpendicular and parallel flow configurations had an inverse relationship with the difference in the percentage of flow contact area for the aperture distributions analyzed. Thus, the lower the difference in JRC values, the more flow contact area expected in the perpendicular flow direction, which will in turn lead to more heat extracted from the rock. From the variogram model parameters, it was deduced that high geometric anisotropy results in high differences in thermal drawdown and consequently high difference in energy extracted, though the values of the geometric anisotropy do not indicate which direction will give a higher thermal drawdown.
From the field scale study investigating aperture distributions with predetermined geometric anisotropic ratios, the perpendicular flow configuration resulted in improved thermal performance over the parallel flow configuration by about 64 % one year into the simulation and 56 % five years into the simulation. However, for the aperture distributions with geometric anisotropic ratios of 2, which had Hurst coefficients similar to values found in real fractures, the perpendicular flow configuration resulted in improved thermal performance over the parallel flow configuration by 70 %.
The results of this study suggest that placing an injector well in the direction perpendicular to shear or slip may result in favorable thermal performance over a parallel flow configuration.