Modeling of multiphase mass and heat transfer in fractured high-enthalpy geothermal systems with advanced discrete fracture methodology

Fracture networks are abundant in many subsurface applications (e.g., geothermal energy, water resources). These networks often have a very complex structure which makes modeling of flow and transport in such networks slow and unstable. Consequently, this limits our ability to perform uncertainty quantification and leads to increased development and environmental risks. This study provides an advanced methodology for simulation based on Discrete Fracture Model (DFM) approach. Changes to the topology of the fracture network reduce computational complexity while preserving the accuracy of the DFM approximation. The preprocessing framework results in a fully conformal, uniformly distributed grid for realistic fracture networks at a required level of precision. The simplified geometry and topology of the resulting network are compared with input (i.e., unchanged) data to evaluate the preprocessing influence. The resulting mesh-related parameters, such as volume distributions and orthogonality of control volume connections, are analyzed. Furthermore, fluid-flow-related changes introduced by preprocessing are evaluated using a high-enthalpy two-phase flow geothermal simulator. The simplified topology directly improves meshing results and, consequently, the accuracy and efficiency of numerical simulation. The main novelty of this work is the fully open-source framework based on graph theory that simplifies the topology of the fracture networks and explicitly resolves the small-angle intersections within the fracture network. Augmenting the framework with a rigorous analysis of changes in the static and dynamic impact of the preprocessing algorithm, we demonstrate that explicit fracture representation can be used together with maintaining computational efficiency enabling their use in large-scale uncertainty quantification studies.

The paper is organized as follows. We start with the description of the input data 151 used in this study followed up by the theory for preprocessing, topology analysis, and 152 fluid flow and energy transport modeling. Next, we describe all essential ingredients of 153 the proposed framework, including intersection, node merging, straightening, and remov- (i.e., the same pattern exists at several length scales) as discussed in  sos (1995). 197 For the dynamic analysis, it is assumed that the two realistic fracture network mod-   The set of all vertices of graph G is denoted as V = V (G), while the set of all edges 221 of the graph G is denoted as E = E(G). Edges of a graph join two vertices i and j such 222 that (i, j) ∈ E(G) and i, j ∈ V (G). If (i, j) ∈ E(G), it implies that i and j are adja-223 cent vertices of G, and i and j are incident with the edge (i, j) .  3. Adjacency matrix: A(G), which is a square n × n matrix, where n is the num-235 ber of vertices of the graph G. As previously mentioned, if the pair of vertices (i, j) ∈ 236 E(G), they are said to be adjacent. Hence A ij = 1 if vertices i and j are on the 237 edge (i, j). Furthermore, for our purposes, it is assumed that the main diagonal 238 is zero (i.e., A ii = 0), which implies that no nodes are connected to itself. Note

244
A typical input data array F that describes the fracture network contains the pair-245 wise x-and y-coordinates of each fracture segment in the network. The first step is to 246 convert this array into two different forms: an array that contains all the unique vertices 247 in the graph (i.e., V ) and the incidence matrix (B). This is done by using Algorithm 1 248 which is found in the Appendix. An important assumption of this construction of V and 249 B is that no subsegments can intersect in other places than the vertices of the partic-250 ular subsegments. As mentioned before, this is often not the case in the manual inter- indicate intersecting fracture segments. This is illustrated in Figure 3. Translating the 265 type of nodes to the graph notation, node I is of degree one, node Y is of degree three, 266 and node X is degree four. In general, it is unusual that more than two lines intersect 267 at exactly one point. However, our preprocessing method merges nodes and causes sev-268 eral nodes to have a degree > 4. This causes us to consider all intersections of node de-269 grees larger than four to be X type of nodes. This is used to plot the results in a ternary 270 diagram (as shown in Figure 3). Classifying the fracture networks topology in this way 271 allows us to use a proxy for the connectivity. Connectivity is often defined in this con-

276
where N I is the total number of I nodes, N Y is the total number of Y nodes, and N X 277 is the total number of X nodes, we use

279
where w i are the weights, N i the total number of i node types in the network, and i =

280
{Y, X, X+, X + + . . .}. X+ represents a vertex which is one degree higher than an X 281 (i.e., degree five instead of four), X + + two degrees higher, etc. The weights are de-282 termined by the number of lines (but not edges) involved in the vertex type (e.g., N Y 283 involves two lines therefore w Y = 2 while N X++ involves six edges and hence three lines 284 therefore w X++ = 3). Please also note that all vertices of degree two are not used nor 285 important in this analysis (i.e., a curved or a straight line are topologically the same).

286
An alternative connectivity measure is obtained by using the Discrete Laplacian 287 of the graph. This matrix can be used for finding spanning trees of a given graph (i.e., 288 connected fracture sets in the fracture network). Notably, each element of the Laplacian's 289 null-space rational basis describes a connected component of the graph (Spielman, 2010).

290
With this basis, we can find the number of connected fracture sets in our network and 291 also each fracture that belongs to these components (i.e., sub-graphs). The connectiv-292 ity measure is then calculated as the ratio between the cumulative length of the fractures 293 in the largest spanning cluster and the cumulative length of all the fractures in the full 294 network.

295
Geometrical properties such as angles and lengths of the fractures are obtained us-296 ing simple trigonometry rules. An easy and fast way to calculate the angles of a frac-297 ture w.r.t. the x-axis is to decompose the fracture into two components (i.e., ∆x = x 2 − 298 x 1 and ∆y = y 2 − y 1 ). Then, the angle can be obtained using the following equation

300
where is a small perturbation to prevent the case of ∆x = 0.

311
These steps are illustrated in Figure 4 and thoroughly explained in the following sections.

312
Lastly, an aperture correction procedure is proposed, which deals with variable apertures 313 and connecting previously disconnected fractures.  the l f , the more precise the preprocessed network represents the raw data. However, small l f means that the subsequent steps in the algorithm take substantially more time. The two edges intersect directly whenever 0 < t, s < 1 is true (note: < instead 339 of ≤ indicates that the intersections at the end-points of segments are excluded). This 340 simplifies to solving a 2 × d system of equations for each intersection, such as The actual point of intersection is calculated by plugging the t that is obtained from The node merging algorithm, in essence, is sequential. Each vertex (i.e., node) is 364 added to the domain that if it doesn't violate the algebraic constraint. This means that 365 the distance between the newly added node and any other node already in the domain 366 must be larger than l f · h, the node is merged into the closest node already in the do-

382
The length of each fracture segment, L ∈ R m , can be calculated in the following 383 way: Then we define the order of adding segments, O segm , from largest to smallest: where a i is the aperture of fracture segment i, and A is the list of all fracture apertures.

388
From now on, for simplicity, it is assumed that h = 1/2. This means that In order to incorporate low connectivity and fracture networks with variable apertures, a fracture aperture correction is applied during the cleaning procedure. Connected fractures are treated analogous to resistors in an electric circuit. Resistance is equal to the inverse of the hydraulic conductivity, which in turn is a function of the square of the fracture aperture. Figure 5 displays the two different corrections (Type 1 vs. Type 2). Type 1 corrections result in the overlap of two segments after merging of vertices. Type 2 is subdivided further into 2a, which results in an edge collapse, and 2b, which connects two previously disconnected edges. The following equation gives the correction for Type 1R where R i is the resistance of edge i defined as where L i is the length and a i is the aperture of the i-th edge respectively, such that the effective aperture of the corrected edge is given bŷ whereL is the length of the new edge.

435
The Type 2a correction is given byR In the case of Type 2b, an effective matrix aperture is obtained by inverting the permeability of parallel plate flow such that Type 2a Type 1 Type 2b where k mat is the matrix permeability at the location of the particular edge. A further addition to the Type 2b correction is added to preserve the characteristics of impermeable fractures/faults and high permeable matrix, given bŷ where L mat is the gap between the vertices and n is determined by fitting a least-squares solution to tracer simulation on two disconnected fractures with different gaps and characteristic cleaning lengths, given by If n ≈ 1 we have the normal harmonic mean, while n → ∞ is the same as not 436 applying any aperture correction. The parameter n effectively limits the aperture penalty 437 when connecting disconnected fractures, as it was observed from simple numerical ex-438 periments that the aperture correction, in some cases, over-penalizes the effective aper-439 ture. Also, note that n is bounded since L mat can never exceed l f h, since otherwise, these 440 vertices would not apply for merging.

441
To deal with vertices that are connected through a path in the neighborhood of the 442 vertices, a sub-graph is extracted around the vertex that is merged. The shortest path 443 is computed using Dijkstra's algorithm implementation described in Csardi et al. (2006).

444
Whenever there is no shortest path (i.e., even in the neighborhood the two vertices re-445 main disconnected), the effective matrix aperture is used for the resistance instead. This 446 is represented in Figure 6 447 Figure 6 illustrates an essential feature of the aperture correction. Merging vertex 448 12 into vertex 11 results in a reduction of the aperture of the edges connecting vertex 449 10 and 12 as well as 12 and 13. This is undesirable because we want to preserve this con-  In order to evaluate the dynamic performance of the preprocessing algorithm, sev-460 eral flow scenarios are considered, for which the governing equations are specified here.

461
The conservation of mass, in general form, is written as where φ represents the porosity, x cp is the molar mass fraction of component c in phase 464 p, ρ p is the density, s p is the saturation, and q p is the source term of the p-th phase re-465 spectively, and v p is the velocity of the p-th phase. The Darcy velocity of the p-th phase 466 is given by where k r,p is the relative permeability, µ p is the viscosity and p p is the pressure of the 469 p-th phase respectively, K is the permeability tensor, and g is the directional gravita- The following equation describes the conservation of energy required for the geother-474 mal simulations: 476 where U p is the internal energy of fluid phase p, U r is the rock internal energy, h p is the 477 enthalpy of phase p, κ is the thermal conduction, and T is the temperature. All govern-478 ing assumptions and properties can be found in (Wang et al., 2020a(Wang et al., , 2021

529
Following the definition of those two parameters, there is a significant distinction 530 between the two preprocessing strategies described below. The first approach is defined 531 as the "clean" strategy. In this approach, the preprocessing algorithm is executed once 532 with a l f = 1, θ a,min = 18 • , and θ s,min = 2.5 • . The l f remains unchanged in the 533 clean strategy for subsequent coarser meshing results. The second strategy is denoted 534 as the "optimal" strategy. In this strategy, for each subsequent coarser model, the pre-535 processing algorithm is executed with l f = l m . This means that the fracture network 536 in the "clean" strategy remains unchanged when coarsening the mesh. In the "optimal" 537 strategy, the fracture network changes when constructing the coarser models.

569
Besides the angle distribution, it is also important to look at connectivity and in 571 particular, the topology changes to the fracture network. Figure 9 shows the topology 572 of the raw and preprocessed fracture networks in the ternary topology diagram (as ex-573 plained in Figure 3). A large deviation between the raw and preprocessed data is observed, approximately 20% I-nodes, 75% Y-nodes, and 5% X-nodes. Furthermore, with increas-577 ing l f , the preprocessed networks increasingly deviates towards a large X-node percent-578 age (from 5% at l f = 1 to almost 70% at l f = 128 for Whitby and from 10% at l f = 579 2 to 70% at l f = 64 for Brejoes).

580
To illuminate the differences in topology between the fine l f = 1 and the raw data, 581 the degree of the raw and cleaned network nodes is shown in Figure 10. Even after zoom-   Figure 9. A large deviation between the raw data and the processed network's topology in both fracture networks is observed. The reason for this is explained in Figure 10. The Brejoes network converges to the raw data for low l f < 1. The jump in the large l f = 128 for the Brejoes case is expected due to the fracture network becoming extremely coarse. Only a few fractures actually remain, meaning the relative proportion of end-nodes greatly increases.  ing results. Therefore, it seems that the volume distribution and quality of the mesh el-597 ements are improved in the preprocessed results. This is quantified in Figure ,   ing method on these two properties, where mesh quality is a proxy for the orthogonal-620 ity of the control volume intersections. Figure 12 shows the volume distribution as a func-621 tion of l f and l m , while Figure 13 shows the distribution of mesh element quality.

622
The volume distribution obtained after meshing the raw fracture network input is  The mesh element quality obtained after meshing with a small l m behaves simi-633 larly for the raw and preprocessed input fracture data, except for a relatively small amount   (Figure 14 and 15). The water saturation field is shown for 651 the Brejoes network (Figure 16), and finally, the temperature at the production well over 652 time (Figure 17).

653
The boundary conditions and modeling parameters can be found in Table 1 and     volumes for larger l m , contributing to small changes across the scales (see Table 3). The 695 difference between the clean and optimal strategy (i.e., l f = l m ) for small l m (≤ 32) 696 in terms of flow-response is negligible; however, the performance of the optimal strat-697 egy is significantly better.

698
A more significant deviation in Brejoes temperature profile for the optimal case is 699 observed. This is in line with the other observations. This pattern is observed in the an-700 gle distribution in the previous section (see Figure 8). Furthermore, Brejoes fracture den-701 sity is larger (i.e., spacing between fractures is shorter), which leads to a more diffused 702 and less complex temperature distribution. The large connectivity also means a shorter 703 and highly conductive path from injector to producer, resulting in an early cold-water 704 breakthrough.  Brejoes, clean Figure 17. The temperature at the production well over time for optimal (left column) and clean (right column) preprocessing strategies for both the Whitby (top row) and Brejoes (bottom row) networks. Substantial deviation for large l f = lm in the optimal strategy was observed. This does not happen in the clean strategy. This is because the fracture network is unchanged while the mesh is coarsened. This also causes the number of control volumes to remain considerable even for large lm thereby reducing the numerical diffusion (see Table 3 and Table 4).  The numerical performance of the two strategies can be found in Table 3 and Ta ison between the two preprocessing strategies for the conventional nonlinear solver.

734
It is observed that the optimal strategy shows a better convergence in both net-735 works. A reduction in nonlinear iterations of roughly 20% for the coarse models in the 736 Whitby simulations is observed. In the Brejoes simulations, this reduction is almost 45%.

737
The total CPU time for the optimal strategy in the Brejoes network increases slightly  allows handling of realistic aperture distributions and low connected fracture networks. 759 We also contribute a comprehensive investigation of the geometry and topology changes 760 as a function of discretization accuracy and its effect on the dynamic reservoir behav-761 ior. Next, we discuss static and dynamic results of our study and give our recommen-762 dations.

769
The large deviation from the raw topology can be explained through several points.

770
Manual interpretation is usually made in some software (e.g., QGIS) or on the image di-

812
The main idea is that adding a fracture to an already connected network is not nec-  Figure 19). Furthermore, the aperture correction allows us to deal with 820 coarsening of sealing fractures that potentially block fluid-flow, as observed in the syn-821 thetic test case in Figure 18.  This study demonstrates a strategy to simplify complex fracture networks in terms 871 of flow response based on a robust preprocessing approach using graph theory. We show 872 that using raw fracture data for topological analysis and dynamic modeling is unwise and 873 that some preprocessing should be applied to investigate the patterns that exist in the if (x i , y i ) / ∈ V then 5: n += 1 7: 8: if (x j , y j ) / ∈ V then 9: idnodes = nonzero(B(·, idsegms(1))) ∪ nonzero(B(·, idsegms(2))) 15: B(·, idsegms(2)) = 0 16: B(idnodes = k, idsegms(1)) = 1 17: 18: //B and X are updated similarly to Algorithm 4 using "mergelistnodes" and "mergelistsegms" for the removed vertices and edges respectively