INVESTIGATION OF THE EFFECT OF THE DISCONTINUITY DIRECTION ON FLUID FLOW IN POROUS ROCK MASSES ON A LARGE-SCALE USING HYBRID FVM-DFN AND STREAMLINE SIMULATION

Understanding the fluid behaviour in rock masses is of great importance in various rock mass-related engineering projects, such as seepage in tunnels, geothermal reservoirs, and hazardous waste disposal. Different approaches have been implemented to study the flow pattern in fractured porous rock masses. Laboratory experiments can provide good information regarding this issue, but high expenses aside, they are time-consuming and suffer the lack of ability to study field scale mediums. Numerical methods are beneficial in simulating such mediums with the Discrete Fracture Network (DFN) method in terms of costs and time as they offer sufficient flexibility and creativity. In this paper, a Matlab code was extended to study the flow regime in a Dual Permeability Media (DPM) with two point sources in the right and left side of the model as an injector and a producer well, respectively. A high permeability discontinuity with different angles was embedded in a very low-permeability limestone matrix. Pressure equations were solved implicitly with a two-point flux approximation scheme of the Finite Volume Method (FVM). Streamlines were traced in the medium and used to analyse the model’s hydraulic behaviour with the aid of Time Of Flight (TOF) for each point. The results show that the FVM-DFN hybrid method can be used as a fast method for fluid flow in DPM with the aid of streamline simulation to study the fluid flow in a large model with discontinuity.


Introduction
Different types of discontinuities, such as bedding planes, joints, fractures and faults are present in the upper layer of the Earth's crust. Such discontinuities form a complex system that makes significant changes, not only in the mechanical behaviour of the rock mass, but in its hydraulic behaviour as well.
Many problems in today's rock engineering demand a better understanding of the flow process in rock masses. Fluid flow can affect both the hydraulic and mechanical properties of a medium. In many studies, a rock mass consists of an impermeable matrix, which cannot conduct fluids, and a system of fractures is the only channel for fluid to flow through. In many cases of rock mechanics problems, the assumption of an impermeable rock matrix will interrupt the precision of the simulation output. Some previous studies have pointed out that the effective porosity and permeability of a matrix, even at very low magnitudes, can change the flow regime and hydraulic behaviour of a system (Odling and Roden, 1997 overall permeability of the medium (Ružičić et al., 2012). It is clear that understanding the fluids and their behaviour in fractured rock masses, particularly on a large scale, is of great importance in various fields of rock engineering.
Analytical and numerical methods have been developed and used extensively to deal with this subject. Although analytical techniques can provide helpful insight into the flow mechanisms, they are restricted to simple geometry and demand a vast knowledge of advanced mathematics for a more complex fractured medium (Namdari et al., 2016;Wang et al., 2019). On the other hand, numerical methods are flexible in terms of implementation and can be utilized in almost any complex geometry (Seales, 2020).
As a consequence of the importance of this issue, hundreds of research papers have addressed this subject through the last century to this present year. Bear et al., (1993) published a book on flow and contaminant transport in fractured rock. Odling and Roden (1997) simulated the contaminant transport in fractured rock with significant matrix permeability with En-Echelon natural fracture system using Finite Difference Method (FDM).

Zhang and Sanderson (2002) modelled a small Dual
Porosity-Dual Permeability Model (DPM) with 3 parallel fractures under stress. They showed that increasing the maximum principal stress can hugely affect the overall image of flow. Trends, prospects, and challenges in quantifying flow and transport through fractured rocks were studied by Neuman (2005) and it was reported that the main challenge is to quantify the exchange of fluid between fracture system and matrix, which is defined through transfer functions. Baghbanan and Jing (2007) investigated the hydraulic properties of fractured rock masses with correlated fracture length and aperture. Their research paper showed that a suitable ellipse of anisotropic permeability could be reached only in cases with relatively large dimensions compared to the fracture medium length.
With significant advances in computing capacities over the years, numerical methods were used frequently by the authors. Namdari et al., (2012) contracted a comparison between a model without matrix permeability and a model with a very low-permeability matrix with simulating over 1800 realization models through the Monte-Carlo simulation technique. It was shown that even in low magnitudes, matrix permeability can change the whole flow regime and overall permeability of a medium. Siavashi et al., (2014) developed a 3D streamline-based simulation of two-phase flow in heterogeneous porous media. The speed of this method was reported to be much higher than other peer methods in heterogeneous media simulation. In 2015, a streamline tracing algorithm in heterogeneous porous rock to predict Time Of Flight (TOF) and breakthrough curves was studied by Nunes et al (2015).
The equivalent Discrete Fracture Network (DFN) for highly fractured rock mass was reported in the study of Ren et al (2017). They simplified the highly fractured rock to a sparser system of fractures, which decreases the amount of time needed to simulate a model. In 2019, Zhang et al investigated the complexities of the hydraulic fracturing process in naturally fractured rock masses. The Synthetic Rock Mass (SRM) method was used and extended through the combination of the Discrete Element Method (DEM) and the Discrete Fracture Network (DFN) to assess the uncertainties. In 2020, Wang et al. performed a tracer test and streamline simulation for geothermal resources in Cuona of Tibet to measure the impact of injection wells in the production of sink wells in different well patterns of this geothermal reservoir. A streamline simulation of experimental studies on waterflooding under given pressure boundaries was conducted as a coupled Hydro-Mechanical (HM) process by Zhang et al. (2021). With the help of the streamline simulation, the capillary effects at the fluid front were quantified as a function of time and coordinates.
One of the key concepts to address here is that discontinuous rock masses do not behave as a continuum (Nam- Hence, the desired method to treat this issue should be fast enough to be able to simulate numerous models in a limited amount of time. In this case, streamline simulation is a scientifically approved method to investigate the hydraulic behaviour of a heterogeneous medium (Thiele, 2002 In this study, the effect of a large discontinuity with different angles in a low permeability limestone matrix is investigated in a 2D square of 100 × 100 m 2 , which is constructed with 2500 cells in the Cartesian coordinate system. To research the discontinuity effect on pressure and flow paths of fluid, the direction angle of the discontinuity is changed from 0 to 11, 45 and 90 degrees with respect to the X-axis. The low permeability rock mass cells are replaced with high permeability cells representing the discontinuity. Through this approach, the main model consists of two main parts: rock matrix cells, which have low porosity and low permeability, and discontinuity cells with a relatively higher porosity and very high permeability. Within this approach, the discontinuity, which is a line, will be replaced by some consecutive porous and high permeability cells in the same direction as the original discontinuity. The newly employed approach in this paper is to use a hybrid approach of the Finite Volume Method (FVM) and the Discrete Fracture Network (DFN), namely FVM-DFN to form a cell-based method. Streamline simulation is used afterwards to trace the path of fluid particles to give a more detailed vision into the flow mechanisms in a fractured porous media, which is lacking in the previous studies, particularly in the rock mechanics research area. Resulting from this, besides the improvement of runtime (which has been decreased to less than one minute for each simulation, in contrast to DFN-DEM approach with hours of simulation time), the amount of time will not increase even for a large number of discontinuities, because the number of equations to solve will not increase, and the cell numbers remain equal even with a high number of discontinuities through this approach. Moreover, this method provides some additional information against DFN methods. A streamline's path and Time Of Flight (TOF) contours will yield a more detailed insight into the hydraulic behaviour of a discrete model, which are missing in DFN methods. In addition, stagnant regions and the overall flow paths will be obtained as opposed to discrete methods.

Methods
In this section, the description of the model is presented. The basic concepts of the finite volume method (FVM) and the Two-Point-Flux-Approximation ap- proach (TPFA) are presented. The streamline tracing method is discussed afterwards.

Model description
In real subsurface conditions, discontinuities tend to have irregular shapes with generally rough surfaces. Also, sedimentary rocks are not homogenous and isotropic due to the layering nature of these rocks and the stress path during the consolidation process. To study the subsurface flow phenomenon, some simplifications may be applied due to the complex behaviour of the ground. Some widely applied simplifications are the assumption of an isotropic rock matrix and representation of discontinuities with a straight line, which is used in DFN approaches .
The models consist of a 2D horizontal 100 × 100 m 2 square with a large discontinuity and impermeable boundary. The rock matrix is a compacted low permeability isotropic limestone. Two source and sink points exist in the middle left and right edge of the square. The discontinuity exists in the square with four different angles of 0, 11, 45, and 90 with respect to the X-axis with the names F0, F11, F45, and F90 (see Figure 1).
In this study, the rock matrix assumed to be a homogeneous and isotropic low permeability limestone. The boundaries cannot permeate the fluid and are defined as impermeable. Hence, the only source and sink points of the model are injector and producer points. Moreover, it is assumed that the discontinuities are a completely straight and homogeneous porous material with very high permeability compared to rock matrix permeability. The matrix's effective porosity is equal to 5% with 0.1 mD (millidarcies) permeability, while the discontinuities are defined as a porous medium with 20% porosity and a permeability of 800 mD. This scheme aims to investigate the effect of different directions of discontinuity with high permeability in a large scale in the absence of gravity and body forces.

Finite Volume Method (FVM) and Two-Point Flux Approximation approach (TPFA)
In the FVM, the basic mathematical idea is to convert volume integrals in a partial differential equation that contains a divergence to surface integrals, using the divergence theorem (Sokolova et al., 2019). Fluxes at the surfaces of each Finite Volume will be derived afterward. These terms are then evaluated as fluxes at the surfaces of each finite volume. Since an incompressible fluid entering a given volume is identical to that leaving the adjacent volume, these methods are conservative. The technique is used in many computational fluid dynamics packages (Moukalled et al., 2016). Finite volume refers to the small volume surrounding each node point on a mesh. Due to the physical concept of FVM, it is a more robust method compared to FEM and FDM in the context of fluid flow.
In FVM, the basic assumption is the physical concept of conservation of mass, or for an incompressible fluid, conservation of volume. In this study, it was assumed that there is a single-phase fluid in the model. Hence, it can be written: (1) Where: v -Velocity vector (m/s) q -Flow rate (m 3 /s) K -Hydraulic conductivity (m/s) p -Pressure head (m) Rewriting Equation 1 will yield the integral form: (2) Where: n -Normal vector on the surface of the control volume To develop a finite volume discretization, rewriting Equation 2 in integral form using a single cell i in the discrete grid and then applying the midpoint rule to calculate the integral will yield: Where is the centroid on . Figure 2 shows a schematic representation of two adjacent cells i and k in a Cartesian system. If it is assumed that the pressure gradient is linear inside each cell, it follows that (Lie, 2019) Where: -Pressure head at the face centroid (m) The transmissibility is introduced here as T parameter. The continuity of fluxes yield that and for pressure . This leads to two equations below: By eliminating the interface pressure πik, it ends up with the following two-point-flux-approximation (TPFA) scheme (6) Where Tik is the transmissibility associated with the connection between the two cells of i and k.

Streamline-based flow simulation
A streamline is a line that is tangential to the instantaneous velocity direction (velocity is a vector, and it has a magnitude and a direction). To visualize this in a flow, one can imagine the motion of a small marked element of fluid. Since the velocity at any point in the flow has a single value (the flow cannot go in more than one direction simultaneously), streamlines cannot cross (Datta-Gupta and King, 2007; Chen et al., 2020). A streamline is the path of fluid particles flowing from the injection well to the production well. It can be defined as the line that is made of points of equal streamline function values. At any time, the tangential direction of each streamline point is identical to the velocity vector.
The most notable study for efficient tracing of the streamlines was performed by Pollock (Thiele and Batycky, 2006). In Pollock's method, the total input and output flux to each boundary is calculated with Darcy's Law (Batista, 2020). With the fluxes being calculated, the scheme focuses on calculating the coordinates of the exit point of the streamline and after that, it computes the amount of time that it takes for a particle to travel along a streamline from the inlet point to the exit (see Figure 3).

(7)
Where: v x 0 -Velocity in the x-direction x = x 0 (m/s) g x -Velocity gradient in the x-direction (1/s) Since , the exit time for each cell can be obtained by integration of the velocity in x-direction (and in the same way for y-direction) for each face for a known entry point (x i , y i ) and exit point (x e , y e ) and some simple algebraic manipulation: It is obvious that the streamline exits from the face with minimum travel time ( ), so the exit point coordinates are calculated by re-solving for x e and y e given the minimum time Δt m : (9)

Results
In this section, simulation results for pressure contours, streamlines, and TOF contours are illustrated. Each simulation takes under 1 minute, which shows a drastic decrease in runtime for a model with dual permeability in both fracture and matrix (Namdari et al., 2012 and 2016).

Contours of pressure
It was noted in the previous sections that the boundaries are impermeable, and the fluid can only be injected or produced in the source and sink points. Pressure contours in the 4 cases are presented in Figure 4.
Pressure contours for the 4 cases are similar in shape and there is no sign of a difference in the shape of pressure distribution caused by different discontinuity angles. There is a significant issue to address here. The pressure distribution may almost be the same in shape, but the pressure magnitude is very different in each case. In the 0 degree angle, the pressure is smaller, practically by 2 orders of magnitude, compared to the three other cases. For the other three cases of F11, F45 and F90, in which the injection is done in the low permeability matrix cell, the pressure range stays almost identical.

Streamlines
In Figure 5, generated streamlines for each case are shown. Discontinuity is drawn with a solid black line. In some cases, such as F0 and F11, streamline concentration on the discontinuity is very high, showing a high permeability effect of the discontinuity. The high concentration of streamlines in these cases represents the high flow region in these two F0 and F11 discontinuities. In the F45 case, streamlines change their direction, align with the discontinuity and then shift back to their original paths. However, the amount of change in this case is smaller in comparison to F0 and F11.
In the F90 case, streamlines don't show any change in the direction and there are no signs of any change in the streamlines path due to the existence of discontinuity.

Time of flight (TOF)
Time of flight contours are illustrated in Figure 6. The direction of the discontinuity has a considerable impact on the TOF contours.
For F0, the fluid flows very fast on the discontinuity line, and it takes the maximum amount of time for the fluid to access other points on the matrix. It can be seen that the high permeability on the line segment between the source and sink is a preferred channel for fluid flow. The TOF contours show that other than cells on the discontinuity line in the F0 model, TOF magnitudes for other cells in the model are about one order of magnitude higher than the same magnitude for F11 and F45 and two orders of magnitudes higher than the F90. Stagnant regions are formed in the upper and lower left of the F0 model, which are presented with a red colour. In F11, the results resemble the F0 model, with high TOF contours in the upper and lower left corners of the model. The smallest stagnant regions belong to the F45 model, with more uniform contours. In F90, the high TOF contours are formed in the right half of the model.
In contour images, the stagnant regions are more important than the absolute magnitudes of TOF for each cell (Lie, 2019). This is because the contours represent the stagnant regions. If the medium is supposed to be saturated with fluid, the injection and production will be on a steady state cycle, so the absolute value for TOF will be in the second order of importance.

Discussion
For F0 (0-degree angle), the pressure magnitude order is two times smaller than the three other cases. However, the pressure for the other 3 cases, 11, 45, and 90 degrees, are in the same order of magnitude. This is because in the case of the 0-degree angle (F0), the injector and producer points are located in the discontinuity line and the fluid (here, water) is injected directly in the high permeability discontinuity. In this situation, fluid can flow directly from the injector to the producer with very low resistance against the fluid particles' movement compared to the matrix. This is the main reason that the overall pressure of the model is smaller by almost two orders of magnitude compared to the three other cases. With a constant rate of injection, the pressure goes up if the permeability of the injecting area is not high enough to permeate the injecting fluid. This shows that the choice of injection spot is of great importance. Also, the pore pressure will change the effective stress and can result in instability issues. For F0 and F90, the pressure contours are completely symmetrical, but for F11 and F45, there is a slight asymmetrical drift in the contours due to the angle of discontinuity.
A streamline's concentration in a given point represents the amount of flow in that region. This is one of the most useful qualitative data points to focus on in simulating a large scale model. When streamlines get close together, it shows a high flow region. In the DEM-DFN method, which is the most accepted method for simulating a model with discontinuity, the flow trajectories do not represent this kind of information (Namdari et al. 2016;. In the F0 case, the streamlines are concentrated on the horizontal line, namely the discontinuity, which shows a high flow region. For F11, the difference in concentration is more pronounced, which shows a high-flow zone on the discontinuity. For the F45 discontinuity, the curve of streamlines changes in the location for F45. After that change, streamlines continue to move on a smooth curve to the producer point. For F90, there is no change in the streamline's path. This can be explained by the fact that the F90 model is completely symmetrical with regard to injection and producer points, so the pressure of all of the points on the discontinuity are equal and there is no hydraulic potential to motivate a flow due to Darcy's law. It can be stated that with an increase in the angle from 0 to 90-degrees of the discontinuity, the impact of the high permeability discontinuity on the streamline's pattern will decrease. For F0 and F11, the impact is very prominent. For F45, the change in the streamlines path decreases, and for F90, there is no sign of a visual difference in the streamline's path. It can be explained by the fact that for lower angles (e.g., F0 and F11), a great length of the discontinuity is in the direction of the flow (which is in the direction of the pressure gradient or a straight line between the source and sink points) and fluid particles tend to reach to discontinuity to travel faster to the sink point. For higher angles (e.g., F45), the angle between the discontinuity and the flow gradient is large, so the discontinuity effects on the streamlines decrease. However, for F90, the discontinuity is perpendicular to the pressure gradient, making it impossible for fluid to flow on this direction considering Darcy's law and the direction of the hydraulic gradient from the source to the sink. In TOF contours, in the F0 case, the fluid flows very fast on the discontinuity line, and it takes the maximum amount of time for the fluid to access other points on the matrix because the high permeability on the line segment between the source and the sink is a preferred channel for fluid to flow. For other points of the model, it takes a long time for the fluid to travel from the injector to the producer compared to the other three models. The amount of TOF magnitudes in other locations rather than discontinuity is between 10-100 times larger compared to TOF magnitudes for the same points in the three other models. It is obvious that the presence of the discontinuity makes this bold difference in TOF magnitudes.
With an increase in the discontinuity angle from 0 to 11 degrees, the patterns of TOF contours do not change very much. Since the fluid tends to flow in the lowest resistance path, TOF contours somehow represent the ease of the flow process at a point. If the TOF parameter for a point is low, it is due to the fact that fluid particles reach this point faster than the other points of the model. are somehow similar, but the TOF magnitudes differ strongly. This kind of information is lacking in DFN-DEM methods as well, since such methods cannot predict the stagnant regions nor the flowing paths of the fluid particle. Also as mentioned before, the DFN-DEM hybrid models are very slow in terms of running speed that can be from hours to days for a simple dual permeability model. In the used approach in this study, the amount of running time is very small, because this proposed method replaces the discontinuity with porous high permeability cells. So numerous discontinuities can be added to the model without any change in running time.

Conclusions
The main challenge in understanding the fluid behaviour in rock masses is to simulate large scale models within a relatively short time in various engineering contexts, like seepage control in underground excavations, geothermal reservoirs and sites of hazardous waste disposal. Different methods have been implemented to investigate the flow pattern in fractured porous rock masses. Many previous studies either consider the fractured medium as an equivalent continuum or simulate the medium with DFN approaches with impermeable rock matrix. In reality, almost every kind of rock has some level of permeability. Using embedded discontinuity in a porous and permeable medium is of great interest in many recent studies. The main problem with these studies is that the simulations need excessive time and computer capacity and each simulation can take hours and even days.
In this study, a novel hybrid FVM-DFN approach is used. The discontinuity is replaced by high permeability cells in different angles of 0, 11, 45, and 90 degrees within a Cartesian network of low permeability limestone matrix of 100 × 100 m 2 . Using a Two-Point-Flux-Approximation (TPFA) method, fluid injected from the left edge of the model (injector or source point) and produced in the opposite side (producer or sink point). Results are presented through pressure contours, streamlines, and Time Of Flight (TOF) contours. The results show that for F0 and F90, the pressure contours are completely symmetrical, but there is a slight asymmetrical drift in the contours due to the discontinuity's angle. Also, for F0, the pressure range in the model is almost two orders of magnitude smaller than the other cases, which can be justified by the fact that injection and production points are on the discontinuity layer and the pressure will be decreased compared to other cases according to Darcy's law. Streamline paths show that for smaller discontinuity angles of F0 and F11, the concentration of streamlines are more prominent due to the fact that the discontinuity is the main channel or flow in these cases. In contrast, for F45 there is a small bending in a streamline's path. For F90, the discontinuity doesn't change the pattern of streamlines, because fluid does not flow through the discontinuity, due to the fact that the pressures of all of the points on the discontinuity are equal and there is no hydraulic potential to motivate a flow due to Darcy's law.
TOF contours illustrate that for F0 case, fluid forms stagnated regions in the upper and lower right corners of the model, because the discontinuity is a highly preferred path for flow and in other points of the model, fluid flows slowly and the stagnation initiate. For F11 and F90, stagnant regions are formed in the model's corners, whilst the stagnant regions in F45 are very small compared to the other cases.
With combining all the contours of pressure, streamlines and TOF, it can be seen that the direction of discontinuity can change the flow regime in DPM rock mass. This is very important in drilling, tunnelling, geothermal reservoir management. These results show that the pressure and flow rate magnitudes are not enough to interpret flow in a rock mass.
From a technical point of view, this study also reveals that runtime has drastically decreased with this hybrid FVM-DFN method in comparison to DFN-DEM approaches. With every run, streamline simulations combined with pressure and TOF contours can provide a detailed insight into the flow pattern of a discrete model with a porous matrix, which can help the engineers in different fields to understand how to operate the projects more efficiently. Also, stagnant regions are hard to detect in DFN-DEM approaches because these approaches break the rock matrix to numerous ultra-small fractures. With these methods and the combination of streamline and TOF contours, the detection of stagnant regions will be straightforward in comparison to other methods.