Hydraulic simulation in fractured media using 1D equivalent pipes-FratCond tool module II

Hydraulic simulation in fractured media is a complex modeling, since it deals with formations that have anisotropic characteristics. Volcanic and metamorphic rocks are examples of materials that make up fractured formations, which commonly have two distinct regions, in terms of fluid movement.. The fractures (or joints), that are planes along which stress has caused partial loss of cohesion in the rock, representing a plane of weakness (discontinuity) in the rock (SINGHAL and GUPTHA, 2010). These planes are quite permeable, facilitating flow occurrence. Meanwhile, the rocky matrix, formed by the cohesive material, presents low or no permeability, hindering the passage of fluids.

Considering the third model category, some studies carried out by different researchers (CACAS et al., 1990a;SEGAN andKA-RASAKI, 1993;DERSHOWITZ, 1996;UBERTOSI et al., 2007;XU et al., 2014) noted the potential of good flow modeling in fractured media, using 1D pipes and channels. This intention aims to reduce computational effort involved in this simulation, considering that such phenomena generally involve complex differential equations, which can be simplified, using flow balance equations in the pipe network, formed from the fractures. In addition, it is known that in fracture planes, only a few regions are effectively wetted, while others remain dry, due to the local variations in aperture and roughness in fracture surface, causing the channeling phenomenon, reported in several studies involving fractured media. According to Ubertosi et al. (2007), in general, only 30% of the fracture surface area is covered by flow, forming preferential paths.
The way to obtain conduits and channels representing discrete fractures is visualized in literature in two modes: from a single representative path, which interconnects the fracture center (considering the different geometric shapes used in its representation) to the intersection central points with the neighboring discontinuities (CACAS et al., 1990a;CACAS et al., 1990b;MORENO e and NERETNIEKS, 1993;DERSHOWITZ, 1996) or by building a conduit mesh in the fracture plane, with variations in the properties that influence fluid transit BODIN et al., 2007). In these studies, both permanent flow simulation and solute transport are presented with satisfactory results, when calibrated and compared with field data.
Considering this context, this study aims to present the aspects involved in the development of a computational tool for permanent flow simulation in fractured media, considering the use of equivalent 1D equivalent pipes. This tool corresponds to FratCond tool module 2, being responsible for the 1D pipe network generation and for the hydraulic heads and flow rates determination, considering contour conditions inserted by the user. The geometry is based on stochastic fracture generation, obtained from module 1 (REIS, ALAMY FILHO, 2018). After the simulation, the tool has graphs and worksheets available for results visualization and analysis.

METHODOLOGY
FratCond module 2 was developed using MATLAB program ming language, being coupled to module 1, as presented by Reis e Alamy Filho (2018). In order to execute module 2, it is necessary to execute the first one, since all model geometry depends on the results obtained in module 1. This generation is based on the main probability density functions (PDFs) defined for the variables that characterize the fractures: length (lognormal distribution), aperture (lognormal distribution), orientation (dip and direction diving; Fisher distribution) and location (fracture center; uniform distribution). From statistical values of these variables (mean and standard deviation), these PDFs can be used to generate fracture networks with similar characteristics to the original, allowing the hydraulic modeling of such media (REIS, ALAMY FILHO, 2018). More details on stochastic approaches related to fractured media can be found in Cacas et al. (1990a), Cacas et al. (1990b), who modeled the Fanay-Augères (France) fractured system using this concept; and also Jing (2003), Dershowitz, Pointe e Doe (2004) e Chilès (2005), where complete reviews about this method can be found.
The module 2 is responsible only for the hydraulic simulation itself.The hydraulic simulation proposed by Fratcond tool is based on the reduction of the discrete fractures network, stochastically generated, to a network of equivalent 1D pipes, as proposed by Cacas et al. (1990a), Cacas et al. (1990b), Dershowitz (1996) and Outters et al. (2000).
These pipes are generated, in parts, from the connection between the intersection midpoint of two fractures and their respective centers. Figure 1 illustrates this generation process, considering fractures represented by Baecher model (circular disks). Each pipe (e.g. 1 , 2 , 3 in Figure 1) is divided in two parts (e.g. 1 and 1 in Figure 1), once each part occurs at the surface of a different fracture, that form the intersection analyzed, with a corresponding change in its orientation, in terms of direction and dip. The repetition of this point interconnection, for all fractures connected to each other, allows obtaining the pipe network representative of the fractured system under study (CACAS et al., 1990a;CACAS et al., 1990b;DERSHOWITZ, 1996;OUTTERS et al., 2000).
It is important to note that the pipe reference nodes are always the fractures center points (black dots in Figure 1), where the representative hydraulic load values will be determined for each fracture. In addition, different pipes will meet in these points, demanding a flow rate balance. The intersection midpoints (red dots in Figure 1) are only basis for pipe network tracing, but do not configure a node itself.
Every fracture that presents intersection with a second one has its representative pipe traced. Isolated fractures, which do not intercept with others, will not participate of pipe network, since there is a dependence on the intersection existence for the pipe tracing.  (2020) For FratCond module 2 execution, two additional information must be provided by the user as input data, which enable hydraulic equations resolution. The first data are the hydraulic head values in the volume simulation limits, which function as a contour condition for the hydraulic simulation. At least one of the four faces that are perpendicular to the xy plane ( Figure  2), forming the cubic or prismatic volume considered in the simulation, need hydraulic head information. In this way, a regional hydraulic gradient can be obtained, indicating flow trend direction inside the simulated formation. Such contour conditions can be acquired in the field, using water level of observation wells, in limits close to the volume considered for the simulation; or even estimated by the user, from observation wells data, that are not necessarily in the perimeters of that volume.  (2020) The second input data required is conductance. Conductance is a coefficient of proportionality between flow and hydraulic gradient, in each one-dimensional equivalent pipe (Moreno and Neretnieks, 1993;Dershowitz, 1996). Two approaches are possible: considering mathematical relationships based on fracture features, such as transmissivity and fracture opening (Dershowitz, 1996); or from the use of statistical values and a lognormal distribution. In the last case, statistical values are obtained from estimation or from simulated model's calibration with field data (considering the relationship between flow, conductance and hydraulic gradients), once conductance is not a physical parameter, directly determined in the fractured media.
Prior to the hydraulic calculation itself, some actions are still necessary for successful simulation. The first action is to separate the equivalent 1D pipes group with more connections, in relation to secondary groups, with less extension or isolated fractures. It is observed in real fractured media that some fractures groups are isolated and do not contribute effectively to water flow. The same happens in the simulations executed in FratCond tool, where it is impossible to consider all the fractures generated in module 1 for hydraulic simulation, which would require multiple contour conditions and could reduce fidelity, when comparing model to reality. Thus, the chosen option was to keep only the longest pipe group, which have a greater chance to traverse the whole simulation volume. This isolation was performed using graph function, available in latest MATLAB versions (later than R2015b). This function analyzes a two columns matrix, considering elements in the same line as nodes. Then, the function informs which paths, formed by these nodes, are connected to each other, through a connection identifier. The identifier with the highest number of repetitions represents the longest pipe group, among all other paths previously obtained.
Another necessary action is to carry out contour conditions, from volume simulation faces to the most extreme pipe network nodes. This transport needs to be done, once the pipe network limits do not always match the volume simulation limits. Thus, the contour conditions in pipe network will be approximated from the values inserted by the user at the volume simulation faces This approximation was made using linear interpolation, considering gradients formed along the x and y axes, which are limited by the volume faces that have the insertion possibility of initial hydraulic heads by the user. Figure 3 illustrates the variables involved in this transport along the x-axis of the simulation volume. Along the y-axis, the process is analogous. The respective meanings of each term are: After completing these actions, the next step is to solve the system of mass balance equations (Equations 1 and 2). Here, we seek to calculate flow rates of each fracture (and in each 1D pipe), in a steady state, as well as total head values in each pipe system nodes, which corresponds to fractures centers. Conductance and pipe lengths are constants and defined previously to this resolution. It is important to know that inflows and outflows are not considered in these equations. where: : flow rate in 1D equivalent pipe, formed between fractures i and j (m³/s); : conductance between fractures i and j (m³/s); : hydraulic gradient between fractures i and j (m/m), given by total head difference (m), divided by pipe length (m).
In addition, the boundary conditions transported to extreme nodes of pipe system must be considered. In this way, total head in the extreme node is equated to the value calculated in transport, as previously described. At least 1 and up to 4 additional equalities can be inserted, depending on the number of boundary conditions inserted in the simulated model. To solve this system, the LSQR method (Least Squares with QR Factorization), developed by Paige and Sounders (1982) and available in MATLAB, was used. A maximum residual error of 10 −4 and a maximum iterations number of 10 6 was considered. Such values were considered sufficient to obtain good results with the use of this numerical method, without the occurrence of numerical instabilities. With the system solved, the values of total head are obtained in each of pipe network nodes. For consequence, we can calculate values of hydraulic gradients and flow rates, replacing such results in the previously presented equations.

RESULTS AND DISCUSSIONS
The results below show the products obtained with the implementation of the method presented previously. FratCond tool user can exploit all these products. All results shown here take into account the fractures generated stochastically by FratCond module 1, from a proposed example, as presented by Reis e Alamy Filho (2018).
As discussed previously, additional input data are required for the proposed hydraulic simulation execution. The user can input this data using the module 2 graphical interface (in Portuguese), illustrated in Figure 4. The blank spaces should be filled with the total head in simulation volume faces and the user need to choose the way conductance will be determined, in the model to be simulated. After this, the simulation can be performed, from the activation of the "Simulação Hidráulica" (Hydraulic Simulation) button. After a few seconds, the user receives a confirmation of successfully routines execution. From this moment, the user can explore simulation results, using the section "Explorar resultados" (Exploring results), which has a list of 6 graphs and 1 spreadsheet, which can be exported. To do this, the user choose the result of interest and then click on the "Ok" button, so that it is displayed. Figure 5 presents the first result available in FratCond module 2, with the 1D pipes network, generated from the fractures obtained stochastically, in the hypothetical example presented by Reis e Alamy Filho (2018). The image shows the fracture identifiers (numbers in red) next to their centers (pink dots), as well as the points of intersection between the fractures (black dots) and their respective identification numbers. The blue lines represent the equivalent 1D pipes. We noted, as expected, the sequence of fracture center-intersection-fracture center interconnection respected, following the methodology used. We can also see that not all fractures are connected to each other, there are pipe groups isolated inside the simulation volume. Secondary pipe networks can be seen near the left side of the volume (between nodes 8 and 25) and the right side (between nodes 15 and 55). The pipe group with the highest number of connections occurs in the lower region of the simulated medium, bounded by nodes 11 and 27. This section will be isolated and hydraulically simulated. Source: Authors (2020) Figure 6 shows the graph available for conductance results analysis, considering all equivalent pipes generated by the model or only the pipes belonging to the longest group. This parameter works as proportionality coefficient between flow and hydraulic gradient. In this hypothetical example, the conductance was defined from fractures properties, according to the methodology of Dershowitz (1996). Results can be checked from the displayed color scale. Light green tones indicate lower values, while dark blue tones indicate higher values. These values are changed at the pipe nodes, indicated by the pink points (fracture centers). The identifiers now shown are pipes identities.
Observing the results obtained in this proposed hypothetical example, we note a high variation of this parameter, when determined from the fracture features. Among them, the one that cause the greatest influence on conductance values are the fracture aperture, indicating more or less ease for flow; and how long is the intersection between the fractures that generate each pipe, which is also indicative of ease for fluid passage.
The following figures show direct results of the hydraulic simulation itself. Figure 7 presents results obtained for hydraulic gradients, while Figure 8 presents flow values obtained for each pipe, allowing a preview of the water volume that transits between the fractures. In both graphs, we can see the highest values occur in pipes 14 and 15, which are highlighted by the blue color. This result is consistent, since both quantities are directly proportional. In addition, for this hypothetical example, boundary conditions entered create a larger gradient in the xaxis direction, when compared with y-axis direction. Considering that pipes 14 and 15 are practically parallel to the axis of higher hydraulic gradient, there will be a tendency of greater flow through these pipes.  (2020) It is noteworthy the possibility of inconsistencies in the displayed graphs, due to the precision entered in the numerical calculation and possible rounding performed, especially regarding the flow balance in each pipe node. As an example, looking at pipe 16 in Figure 8, despite its continuity with pipes 14 and 15, this tube had a slightly different flow rate, changing its color class in the graph displayed. However, by looking at the result spreadsheet with the numerical values, it is possible to verify properly that this is an acceptable rounding problem, at the fourth decimal place. Source: Authors (2020) Finally, the last available graph (Figure 9) shows the total head values obtained in each of the pipe network nodes (fracture centers). This graph shows color-scale circles that indicate the respective total head results. Identifiers for fracture centers (red numbers) and pipe stretches (black digits) are also dis-played in this view, and the numbers surrounded by a rectangle indicate pipe nodes that have received the boundary conditions, i.e., they are extreme nodes, closest to simulation volume faces. For this hypothetical example, as boundary conditions were imposed in all possible simulation faces, four nodes are highlighted, with nodes 27 and 11 receiving right and left faces conditions, while nodes 20 and 14 receive from the bottom and front faces conditions, respectively.
We also note that the total head values obtained respect flow tendency imposed by the boundary conditions initially inserted in the model. Thus, there is a reduction in total head from right to left and bottom to front, with the largest contributions, in terms of flow rates, appearing in pipes that are parallel to the x-axis. Source: Authors (2020) As a last form of results analysis, the tool has an exportable spreadsheet available (Figure 10), which gathers the main data obtained in the hydraulic simulation, in numerical form.
The first tab presents data used to determine conductance, such as pipe length and width, the average pipe aperture, the average transmissivity and the respective conductance value obtained. For a complete understanding of the relationship between these parameters, it is recommended to consult the study of Dershowitz (1996). A pipe identifier, followed by the indication of the respective fractures that form it, orders all these results. The second tab presents the results of the hydraulic simulation itself, for the longest pipe group. At the top, data concerning the entered and transported boundary conditions, the numerical resolution conditions of the proposed system of equations are displayed. At the bottom, for each pipe, it is listed the following data: pipe length, conductance, total head at the beginning and end of the pipe, the hydraulic gradient and the respective flow rate.

CONCLUSIONS
Considering the results presented, we conclude that FratCond module 2 is able to take advantage of the results provided by module 1 and to simulate a fractured medium, considering the approach of equivalent one-dimensional pipes. The process to obtain pipes geometry and the calculation of hydraulic variables were well implemented, considering the methodology used in this study. The results analysis is facilitated using graphs and spreadsheets available to the user. We expect, in a possible study continuation, an additional module implementation, capable of evaluating the solutes transport in fractured media, considering 1D pipes approach; and possible simulations of fractured media with real data, for a full evaluation of FratCond tool potentialities.