Modelling the Flow in the Utah FORGE Wells Disrete Fracture Network

Abstract

The focus of this paper is the efficient numerical solution of the fluid flow in the Utah Frontier Observatory for Research in Geothermal Energy (FORGE) reservoir. In this study, the public data available for Discrete Fracture Networks (DFN) around well 58-32 is used to represent the DFN. In this research, a novel computationally efficient method called Hele-Shaw (HS) approximation is used for modeling fluid flow in FORGE well. An analysis of the influence of fracture intensity in a network is carried out using the HS method. The HS method was validated by solving the full Navier–Stokes equations (NSE) for a network of eight fractures. A good agreement was observed between the evaluated results (average deviation of 0.76%).

1. Introduction

The term Enhanced Geothermal Systems (EGS) describes systems that utilize geothermal resources to produce economical amounts of heating from low permeability and porosity reservoirs. For a long time, experts have seen EGS as a possibility for developing geothermal systems outside of traditional hydrothermal areas [1]. EGS heat reservoirs consist of low permeable, conduction-dominant fractured rock systems [2]. After injecting the fluid, hydraulic stimulation creates a network of fractures with more open connections, which allows hot rocks to provide heat to the fluid [3]. The interaction of fluids with hot rocks can alter fracture permeability or cause new fractures to form. This fluid is then pumped to the surface using a production well.

The US Department of Energy (DOE) has selected FORGE as a test site for scientists and engineers. In order to commercialize energy production from low permeability hot formations, or EGS, this project aims to develop new technologies, techniques, and procedures. Milford, Utah is home to FORGE’s laboratory site, just west of Blundell Power Plant (see Figure 1a,b). The geothermal potential of the site has led many previous researchers to develop conceptual hydrological and geological models of the site. Roosevelt hot springs confirm a source of heat under the mineral mountains.

The FORGE reference DFN was generated by the Utah team using the FracMan 7.7 software [5]. Data from well logs is incorporated into the DFN in order to create a single hydrologic and mechanical system [6]. Fracture-set orientations and intensities were determined from well 58-32 (see Figure 2), whereas fracture sizes were determined from nearby trace length data collected in the mineral mountains. Micro resistivity images and dip data are provided in real-time by the Formation Micro Imager (FMI) in water-based mud that is used to evaluate fracture orientations (detailed information is provided in [7,8,9]). Afterward, the Utah team imported the DFN into Schlumberger’s software program known as Kinetix, which investigates hydraulic–natural fracture interaction. To analyze the heat transfer, the Utah team used a software application called Transport of Unsaturated Groundwater and Heat (TOUGH2) [10].

Hydraulic fracturing involves injecting fluid at high pressure into subterranean rock formations. Research has tended to focus on traditional modeling methods rather than computationally efficient approaches for simulating this process. The FORGE well has been analyzed in several investigations [7,11]. In particular, Finnila et al. [7] focused on the evaluation of fracture intensity, flow rate, and permeability. For this assessment, they have used the FMI method alongside TOUGH2 and Kinetix software. Riahi et al. [12] developed a coupled hydro-mechanical numerical simulation to design a conceptual FORGE well geometry. Simulating fluid flow was accomplished through the use of Cubic Law (CL) method. They draw our attention to the limitations of current modelling methods that are (1) disregarding fracture variations and orientations, and (2) neglecting small fractures. Simplifications are made to reduce the computational costs that result in inaccuracies in the approximation. Nadimi et al. [8] estimated the performance of the reservoirs by modeling the fluid/heat system using TOUGH2. According to the authors, the CPU constraint was the main downside of their approach. Lu and Ghassemi [13] simulated the first experiment of EGS Collab performed by [14,15]. The EGS Collab is intended to establish validations for subsurface models and research against controlled, small-scale in situ experiments on rock fracture behavior and permeability enhancement in the field on a small scale [16]. Their fracture permeability solution is based on the CL. Combining Darcy and CL enabled them to investigate how fracture deformation affects reservoir performance.

Figure 2. Geological map of the Utah FORGE site based on previous works [17].

Figure 2. Geological map of the Utah FORGE site based on previous works [17].

A flowback event is one of the final stages of hydraulic fracturing. As the name suggests, flowback refers to a process in which the fluid used for hydraulically fracturing a shale formation is withdrawn from the well at the surface [18,19,20]. In the drilling process, this procedure is performed in preparation for a subsequent phase of treatment, or to clean the well and transition it to the production phase. The flowback consists of water, dirt, sand, and chemicals. To protect surface and groundwater from contamination, sudden pressure changes must be avoided, and reclaimed fluid must be properly disposed of. Xing et al. [21] numerically simulated the flowback event. They presented three sets of natural fractures, which included a total of 2000 fractures. The simulation method is called the Distinct Element Method (DEM) that is based on the CL. The modelling was conducted on several stages for simulating different zones within the network. For understanding the influence of pressure change on fracture closure, their numerical model analyzed the reservoir pressure change under flowback operation. Note that their model was calibrated based on experimental data. Finally, they demonstrated the impact of rebound pressure continuity on the prevention of fracture closure.

Computationally efficient approaches for modelling fluid flow in EGSs are of interest. Many researchers have attempted to develop computationally efficient approaches. For instance, De Dreuzy et al. [22] simulated single-phase fluid flow in eighteen fractures with ten random configurations. The authors utilized Poiseuille’s law for approximating the network’s head. A group of researchers [23] developed an optimization formulation for simulating fluid flow in DFNs. Since their method was not developed based on the discretization of partial differential equations, the resulting approach was computationally efficient, compatible with using parallel computers and consistent with conventional numerical approaches. This method is a Graphics Processing Unit (GPU)-based approach that is applicable to parallelized computers.

Despite the fact that the DFN data are available to the public, only few researchers have attempted to analyze the fluid flow behavior in the network. The authors have mostly relied on simplified methods, such as CL. The CL is an empirical relationship derived from the NSE for laminar flow in a fracture [24]. The HS approximation is used to describe flow in thin, planar fractures or gaps where the width (aperture) is much smaller than the other dimensions [25]. Although the CL approach is interesting, it suffers from inaccuracy due to overestimation of flow rate and failure to consider the effect of fracture surface roughness [26].

In this paper, the introduced HS approach in our previous work [25,27,28] is employed to analyze fluid flow behavior under different geological conditions in the FORGE reservoir. The Leapfrog software is used to produce high-quality presentations of field data [29]. The COMSOL [30] Multiphysics package v5.6 was utilized to perform numerical simulations. In order to verify the accuracy of the HS method, the HS solution was compared with the NSE solution for a network of ten circular fractures (see Section 3.2). Using the HS method, the influence of the number of active fractures, network intensity, and aperture size were analyzed on the fluid flow behaviour.

2. Methodology

2.1. Numerical Methodology

The solution of the full NSE for a large number of fractures is computationally prohibitive for realistic networks. As an alternative to the NSE approach, the HS method has a well-behaved governing equation that does not require supercomputers for modeling fluid flow in a large-scale DFNs.

The governing equations for an incompressible, Newtonian, single-phase, steady-state flow are the NSE. Equations (1) and (2) represent a set of nonlinear partial differential equations for the velocity and pressure field [31]. These equations can be expressed as momentum and mass conservation, respectively. Accordingly,

$${\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+(\mathbf{u}.\nabla )\mathbf{u}={\displaystyle \frac{-1}{\rho }}\nabla P+{\displaystyle \frac{\mu }{\rho }}{\nabla }^2\mathbf{u}$$

(1)

$$\nabla ·\mathbf{u}=0$$

(2)

In these equations, $\rho$ is the density (kg/m³), $\mathbf{u}=[u,v,w]$ is the velocity vector (m/s), P is the pressure (Pa), and $\mu$ is the viscosity (Pa.s). To fix ideas and because of its relevance in ground water transport, throughout this work we use water at 200 °C as the working fluid at the target depth of the well 58-32 with the following properties $\rho$ = 963.33 (kg/m³), $\mu =1.3\times {10}^{-4}$ (Pa.s).

An approximate version of the NSE can be obtained by considering the flow in two directions and integrating the NSE across the fracture aperture. As the inertia forces are considerably smaller than pressure and viscous forces when the Reynolds number of the flow-through fractures is small, they are neglected. The Reynolds number in fractures can be calculated using Equation (A21) in the Appendix A. Accordingly, the velocity is averaged across the fracture aperture and expressed as the depth-averaged velocity $\mathbf{V}=[\overline{u},\overline{v}]$. Where $\overline{u}$ and $\overline{v}$ represent flow parallel to the fracture wall. This simplification is possible due to the smaller aperture compared to the fracture extent (orders of magnitude). The depth averaging approach reduces the computational costs considerably.

This depth-averaging approach reduces the computational complexity by transforming the domain from three-dimensional with a three-component velocity field to two-dimensional. Consequently, the simplified governing equations for flow between parallel plates separated by an aperture (h) are expressed by Equations (3) and (4). Please refer to Appendix A for the complete set of governing equations.

$$\mathbf{V}=-{\displaystyle \frac{h^2}{12\mu }}\nabla P$$

(3)

$$\nabla .\left(h\mathbf{V}\right)=0$$

(4)

The pressure distribution is governed by the Laplace equation, derived by substituting Equation (3) into Equation (4). Consequently, the HS equation (Equation (4)) can be solved over the fracture’s center plane, significantly reducing computational time.

We solved the governing equations using COMSOL v5.6 [30]. Our approach is similar to that described in our previous work [25,27,28], where we solved the incompressible NSE in a three-dimensional computational domain using the single-phase laminar flow (SPF) module. We applied the HS approximation to the corresponding two-dimensional surfaces using the Coefficient Form Boundary PDE (CB). The steady-state solver was used in order to solve both equations. The NSE and HS were solved using the Generalized Minimal Residual (GMRES) and MUltifrontal Massively Parallel Sparse direct Solver (MUMPS) algorithms, respectively.

2.2. Direction and Angle of Fracture

Understanding fracture features requires the observation of fracture geometry from boreholes or exposed rock faces. DFN fracture features include different orientations and angles. In mathematics, lines and planes in space are defined by equations. However, geophysicists use a different approach for characterizing fracture features. A compass is the main tool for this characterization. There are two different methods for representing the compass direction, azimuth and quadrant (see Figure 3). Note that the azimuth direction is more commonly used in geology.

To characterize fracture features, the International Society of Rock mechanics (ISPM) proposed several parameters that includes strike, dip, trend and plunge. Strike and dip are used by geophysicists to identify planes’ orientation. Additionally, lines orientations are described by trend, plunge and pitch. Below is a brief description of these parameters:

• The strike represents the orientation of a planar feature. The strike can be described as the direction of the red line in Figure 4 as occurred at the intersection of horizontal plane and the plane of interest (the blue plane). The strike direction can be measured by a compass.
• The term dip refers to the angle of inclination of a planar feature as measured from a horizontal datum (see Figure 4). This measurement refers to the angle (inclination) of a planar feature. In other words, the dip angle is the angle between the horizontal plane and the plane of interest. The dip direction can be determined by moving 90 degrees clockwise from the strike direction. The strike direction in Figure 4 is 000, therefore the dip direction will be due east or 90°.
• The trend of the line is a horizontal projection of the line of interest (red line in Figure 5a). The top view of the box in Figure 5a allows us to see only the trend of the line.
• The pitch of the line refers to the angle between the horizontal line and the line of interest (the yellow line in Figure 5b) within the plane.

3. The Use of Hele-Shaw Approximation in Practice

There are usually four steps involved in the construction of a DFN model in a new drilling location [7]. Defining the boundaries of the modeling region is the first step in the modeling process. Furthermore, the fracture network can be constructed within the defined boundaries by using the random generation method described in [27]. Next, changes must be made to the location and orientation of the fracture. Lastly, it is necessary to calibrate the fracture aperture sizes.

To accommodate a geothermal reservoir previous studies have assumed a volume between 1 and 5 km³ [33,34,35] at a depth of approximately 500–4000 m below the surface [36]. In Figure 6, the DFN model region of FORGE is presented as an example. The lithology in that area is classified into two broadly defined units consisting of granitic basement rocks and overlying sedimentary deposits [37].

The FMI log interpolation can be used to adjust the frequency of the fractures in different regions. The FMI is a highly advanced dipmeter tool capable of producing high-resolution resistivity images of borehole walls [38]. With the FMI data, it is possible to perform adjustment on the network frequency and the aperture size. The random distribution method presented in this paper is the main tool for this adjustment. Detailed description of the FMI approach can be found in [7]. Having the configuration of the fractures, the HS approximation can be utilized for simulating the fluid flow and obtaining the hydraulic resistance.

3.1. Mesh and Boundary Conditions

The boundary conditions are applied to the edges resulting from the interaction of fractures with production and injection wells. Constant pressure Dirichlet boundary conditions for inlets and outlets are prescribed. Figure 7a represents the imported DFN geometry around well 58-32 into COMSOL software from the public data provided by the FORGE team [39]. This data contains individual fracture orientations as dip, strike, trend and plunge. This network is in the shape of a cube around well 58-32, with dimensions of $1\times 1\times 1$ km. Several wells are located near well 58-32. The two closest are wells 78-32 and 68-32, which have been specifically drilled for monitoring purposes. To produce geothermal energy, two wells are required: one for injection and one for production. An approximate diameter of 0.3 m is required for the wells [40]. Well 58-32 has a diameter of 0.25 m [41]. The data for the production well was not made available by the FORGE team. Therefore, we have assumed two wells on both sides of the network with a diameter of 0.3 m. Most simulations for this network were conducted with a global pressure difference $\Delta P$ of 0.01 (Pa). All fractures were assumed to have a constant aperture size of 0.01 m in all analyses except for the analysis of fracture aperture size.

To solve the NSE for the verification case, the fracture domain was discretized using a mixture of tetrahedral and triangular mesh elements, for a total number of 68,853 elements. The corresponding two-dimensional HS domain was discretized using 16,868 triangular elements. The mesh network on 350 fractures around well 58-32 is depicted in Figure 8a. The number of elements in this figure is 490,980. To represent the mesh quality, skewness is reflected over fractures (see Figure 8b). The green color of all elements in Figure 8b indicates a high-quality mesh. An average skewness of 0.807 is measured from the mesh. To ensure numerical stability and mesh independence, a mesh sensitivity analysis was performed. The findings of this analysis for all models are not presented here for the sake of conciseness. The flow rate was measured from the cross-section of the network shown in Figure 7b. As presented in Figure 9, the result did not change after the number of elements exceeds 490,980. Therefore, the main analysis of this model was performed with 490,980 elements. The mesh sensitivity analysis demonstrated that the solution converged as the mesh underwent refinement.

3.2. Verification

NSE have been used for validation since they are the gold standard for assessing simplified models [42,43]. Since the verification of 350 fractures by the NSE is computationally prohibitive, we examine the validity of the HS approximation by modeling fluid flow in eight fractures (see Figure 10). In this simulation, Fracture 1 and Fracture 8 are square shapes. The reason for this is the simplification of the meshing process. During this verification, constant geometrical properties were considered, including the aperture size (0.2 m), intersection length (2 m), and angle (0°). The aperture of 0.2 m might appear unrealistic. Constant pressure Dirichlet boundary conditions were prescribed for the inlet and outlet. In this simulation, a global pressure difference $\Delta P$ of 1.0 $\times {10}^{-7}$ (Pa) was assumed. On this basis, the full NSE were solved first. Next, the fluid flow was modeled in the corresponding two-dimensional geometry using the HS approach. Then the pressure distributions obtained by both methods were compared. As presented in Figure 11, the predicted pressure distribution by the HS method (Figure 11a) and the NSE method (Figure 11b) are identical. Figure 11c represents the deviation between evaluated pressure by both methods over Fracture 4. The average difference is less than 0.76%. Accordingly, the satisfactory agreement between the NSE and HS results confirms the validity of this method for modeling the FORGE problem. It is worth pointing out that the computation time required for solving the full NSE for this particular geometry was 12 min, while the solution time for the HS approximation took only 2 s.

It should be noted the deviation of the HS approximation from the NSE was calculated using Equation (A23).

4. Results and Discussion

4.1. Impact of Number of Fractures

The impact of randomly removing fractures from different regions of the network was evaluated. Approximately 20 to 70 fractures were removed randomly, and the flow rate was measured at the cut plane located in the middle of the network (Figure 7). The removed fractures in three realizations are presented in Figure 12 (blue fractures).

Table 1. The effect of removing fractures randomly from different area of the network.

P32	N	QHS (m3/s)	Re	HR (Pa.s−1.m−3)
0.0221	350	4.21 $\times {10}^{-7}$	1.62 $\times {10}^{-3}$	2.374 $\times {10}^{-1}$
0.0194	307	4.14 $\times {10}^{-7}$	1.59 $\times {10}^{-3}$	2.417 $\times {10}^{-1}$
0.0185	293	3.03 $\times {10}^{-7}$	1.16 $\times {10}^{-3}$	3.302 $\times {10}^{-1}$
0.0163	278	1.98 $\times {10}^{-7}$	7.63 $\times {10}^{-4}$	5.043 $\times {10}^{-1}$
0.0148	253	1.73 $\times {10}^{-7}$	6.66 $\times {10}^{-4}$	5.773 $\times {10}^{-1}$

details the evaluated flow rate and hydraulic resistance for networks with different numbers of fractures. For calculations regarding hydraulic resistance, please refer to the Appendix A, Equation (A22). The number of fractures is influential on the hydraulic resistance and flow rate. As anticipated, a larger number of fractures has resulted in a higher flow rate and lower hydraulic resistance. Furthermore, there is a direct relationship between flow rate and network intensity. The more intensive the network, the lower the hydraulic resistance and the higher the flow rate. The intensity of the network P₃₂ can be calculated by dividing the total fracture area by the rock volume. Note that P₃₂ in this analysis ranged from 0.221 to 0.0148. Results obtained for each simulation using the HS approximation have taken 47 s. This fast computational simulation method has enabled us to analyze a real-life problem.

4.2. The Effect of Fracture Aperture

To determine the effect of fracture aperture on hydraulic resistance, the aperture size was varied between 0.001 and 0.01 m. In order to maintain accurate results, a constant global pressure difference of 0.01 (Pa) was prescribed. Note that, the HS approximation provides results with high accuracy at Reynolds numbers below 10 as examined in previous work [25]. The flow rate was measured at the cut plane specified in Figure 7b. The results for aperture sizes ranging from 0.001 to 0.010 are summarized in

Table 2. The impact of different constant fracture aperture on the hydraulic resistance.

h	QHS (m3/s)	Re	HR (Pa.s−1.m−3)
0.001	7.2 $\times {10}^{-10}$	2.8 $\times {10}^{-6}$	138.8
0.002	5.8 $\times {10}^{-9}$	2.2 $\times {10}^{-5}$	17.35
0.003	1.9 $\times {10}^{-8}$	7.5 $\times {10}^{-5}$	5.142
0.004	4.6 $\times {10}^{-8}$	1.8 $\times {10}^{-4}$	2.169
0.005	9.0 $\times {10}^{-8}$	3.5 $\times {10}^{-4}$	1.111
0.006	1.6 $\times {10}^{-7}$	6.0 $\times {10}^{-4}$	0.643
0.007	2.5 $\times {10}^{-7}$	9.5 $\times {10}^{-4}$	0.405
0.008	3.7 $\times {10}^{-7}$	1.4 $\times {10}^{-3}$	0.271
0.009	5.3 $\times {10}^{-7}$	2.0 $\times {10}^{-3}$	0.190
0.010	7.2 $\times {10}^{-7}$	2.8 $\times {10}^{-3}$	0.139

. Due to the larger fracture aperture, the global pressure difference acts over a larger cross-section, thus decreasing flow resistance, resulting in smaller values of corresponding hydraulic resistance.

4.3. The Effect of Fracture Frequency

To examine the effect of fracture frequency on fluid flow in different regions of the network, this study measured the average velocity in eleven cross-sections along the length of the network. In Figure 13, we show a cross-section of the network. To calculate the fracture frequency (P₁₀), fractures are counted at eleven cross-sections and are divided by the length of the network, which is 1000 m.

According to Figure 14, the velocity is consistent with the value of P₁₀ in most areas. Therefore, higher velocities have been observed in regions with lower frequencies. Nevertheless, there are expectations that the velocity should not decrease in response to the network frequency. This is attributed to the fracture disconnections or the occurrence of dead ends. A closer look at Figure 14 shows that cross-section 540 has a higher velocity compared to the fracture frequency. Even though the network is more intense in that area, the velocity has increased.

5. Conclusions

This work has utilized the HS approximation for modelling fluid flow in the DFN around FORGE well 58-32. We compared the results of the HS approximation with the full NSE. The obtained results were in good agreement (average deviation of 0.76% for the flow regime investigated). Using the HS approximation, the computational time was reduced while the accuracy of the result was maintained. For the verification case, we solved the full NSE in 12 min, while for the corresponding two-dimensional geometry, we solved the HS approximation in two seconds. Fluid flow was modeled in 350 fractures using the HS approximation in 47 s. This computational efficiency is the main contribution of this research that has enabled us to analyze a real-life problem. This computationally efficient method for modeling fluid flow in fractures can be integrated with other more advanced models to capture the heat and mass transfer process in fractures such as the one proposed in Wang et al. [44] to achieve comprehensive fluid flow modeling within a reservoir. We examined the impact of the number of fractures and network intensity by randomly removing fractures from different regions of the network. Increasing the number of fractures resulted in less intense networks and lower flow rates. Additionally, the hydraulic resistance decreased with increased network density. The study also examined the effect of aperture on the flow. Based on the evaluated results, the hydraulic resistance decreases as the aperture size increases. The results indicate that fracture frequency directly influences velocity. There is, however, a possibility that in some regions of the network with high frequency, fracture disconnection or dead-ends result in a higher velocity than in other regions of the network with a lower frequency. We believe that the HS method can make significant contributions to the field in several ways. Firstly, it can be coupled with the Darcy equation to model permeability in the surrounding matrix of the fracture network. This coupling does not require massively parallel computers to solve large-scale problems. For large DFNs, using Darcy’s law coupled with the NSE is impractical. The HS method can also be coupled with other equations to model heat transfer or chemical reactions. Further research is needed to determine the influence of the surrounding matrix of large fracture networks on fluid flow and hydraulic resistance. Much remains to be conducted to better quantify the flow in a realistically large and geometry compliant DFN in particular related to the effects of surface roughness, fracture number density and orientation, the coupling between the flow in the fractures and the permeability of the surrounding medium, or the effect of multiphase flow which is characteristic of geothermal reservoirs. These important questions will be the focus of future work.