Numerical fluid flow modelling in multiple fractured porous reservoirs

This paper compares the fluid flow phenomena occurring within a fractured reservoir for three different fracture models using computational fluid dynamics. The effect of the fracture-matrix interface condition is studied on the pressure and velocity distribution. The fracture models were compared based on the variation in pressure and permeability conditions. The model was developed for isotropic and anisotropic permeability conditions. The results suggest that the fracture aperture can have a drastic effect on fluid flow. The porous fracture-matrix interface condition produces more realistic transport of fluids. By increasing the permeability in the isotropic porous matrix, the pressure drop was significantly higher in both the fracture and reservoir region. Under anisotropic conditions in the 3D fractured reservoir, the effect of the higher longitudinal permeability was found to lower the pressure in the fractured reservoir. Depending on the properties of the fractured reservoir, this study can enhance the understanding of fracture-matrix fluid interaction and provide a method for production optimisation.

well as investigating the flow relationship between multiple reservoir fractures. Therefore, this paper aims to develop the reservoir model to include both intersecting and multiple parallel fractures, to produce a more accurate representation of natural reservoir formations. This paper uses computational fluid dynamics to solve Navier-Stokes equations and model the fluid flow and transport in multiple 3D porous fractured reservoirs, under varying input parameters namely pressure, isotropic and anisotropic permeability, fracture-reservoir interaction, fracture geometries, and inertial resistance. The role of the change in fracturematrix interface condition on pressure drop and the flow rate has been investigated in detail.

Methodology
The computational fluid dynamics of single-phase fluid flow and transport of water through fractured porous media was studied using the finite volume method in the commercial software ANSYS FLUENT. The fluid modelling process can be summarised in the following steps: a) Three different fracture geometries were created that would reside within a porous sandstone reservoir. b) An accurate grid-like mesh within the fracture and porous media was generated. c) An accurate model with appropriate boundary conditions was set-up. The three-dimensional fracture-reservoir geometries were developed and simulated for a single phase, steady-state water flow driven by a differential pressure head. The flow was simulated for isotropic and anisotropic reservoir using both wall and porous fracture-matrix interface conditions.

Governing equations 2.1.1 Continuity equation
In an isothermal system the continuity equation for a steady state, incompressible condition can be defined as: (1) where v i is the velocity vector.

Momentum equation
The Navier-Stokes equation was used to model the momentum change in porous media defined in Eq. (2). The Eqs. (1) and (2) are based on isothermal, steady state, incompressible condition assumptions and thus the transient terms are neglected.
(2) where µ is the fluid viscosity, ρ is the fluid density, P is the static pressure, and F i is the source term to account for the flow through porous media, and can be calculated by rearranging the Darcy's Law.
(3) where k is the permeability of the reservoir. For isotropic reservoir, k is assumed to be homogenous, and for anisotropic reservoirs, the effect of change in permeability was studied.

Reservoir fracture geometry
The geometry in this paper was originally developed from an experiment by Karpyn et al. [Karpyn, Alakmi, Parada et al. (2003)]. A core sample of Berea sandstone was artificially fractured using the Brazilian method and scanned using micro-computerised tomography (MCT). This was converted into data points and used to generate a complex digital geometry. Crandall et al. [Crandall, Ahmadi and Smith (2010)] and Nazridoust et al. [Nazridoust, Ahmadi and Smith (2006)] simplified the geometry using the parallel plate method for their studies, which considered water flow through 2D fractures. Fracture images from these papers were used in plot digitiser to generate the geometry using Solidworks. The length of the fracture used in the 3D analysis is 0.0066925 m, as shown in Fig. 1. The reservoir dimensions used for all models were 6.6925×2.5×2.2308 mm (L×B×H). The fracture shown in Fig. 1 was adapted to be slightly more complex, with the addition of arms, as shown in  In the parallel plate model, the fracture is represented by a set of parallel plate openings with different apertures and modelled using Eq. (4) as can be seen in Fig. 2. It is assumed that the total pressure drop in the fracture can be calculated as a sum of the pressure drop in different parallel channels. In addition, the computed pressured drop has to be corrected for the tortuosity. The detail description of the parallel plate model can be found in Nazridoust et al. [Nazridoust, Ahmadi and Smith (2006)]. The pressure drop can then be used to calculate the effective friction coefficient (f), defined as: where ρ is the fluid density, L is the fracture length, Q is the flow rate, θ is the tortuosity, ∆P is the pressure drop, H � is the effective fracture aperture and can be calculated from H avg is the average fracture aperture, and σ is the standard deviation of fracture apertures.
ii) Coarse model As the parallel plate model is a simplification of the complex fracture geometry and in order to capture flow characteristics missed by the parallel plate assumption, highlighted by Petchsingto et al. [Petchsingto and Karpyn (2009)], it is necessary to create a reservoir flow model that is more accurate and representative to real fracture geometries. The coarse geometry was therefore created to represent a real fracture while maintaining the same characteristics (average aperture height and tortuosity) as the parallel plate model. Firstly, the parallel plate aperture distribution was determined from 62 positions along the fracture length. The resulting normal distribution curve is shown in Fig. 3. For each coarse aperture, a value was selected from the curve on the left-hand side (LHS) of the vertical axis and coupled with the value in the mirrored position on the right-hand side (RHS). This process was repeated, and all values were assigned locations across the overall length (6.6925 mm) to ensure the aperture distribution remained constant across the parallel to coarse conversion.
Furthermore, tortuosity has been defined by Crandall et al. [Crandall, Ahmadi and Smith (2010)] as the percentage difference between the matrix and total passage length due to wall roughness and asperities. In order to calculate this value, the passage length had to be determined. Construction lines were positioned across the aperture heights along the length of the parallel geometry; the centre position of each line was connected, and the total length measured. Tortuosity was then calculated using Eq. (6).
where: θ is tortuosity; L e is equivalent length (passage length); L is matrix length.
To ensure the tortuosity remained consistent throughout the coarse fracture the centreline passage was used as a template when sketching the profile of the coarse geometry. This allowed for tortuosity control while using the aperture distribution data. Once the model was generated the tortuosity was checked using the same process as mentioned before, and the overall volume of both parallel and coarse geometry was compared to ensure the same surface area was present. This method resulted in the same fracture characteristics being present in both and meant that the effect of coarse geometry could be accurately analysed. The coarse fracture model is shown in Fig. 4. iii) Z-Varied (using parallel plate model) Natural fractures are very complex due to a number of forces and movement that affect their geometry. Previously only fractures with profiles varying in the X and Y directions have been considered. In this section, the fracture aperture has been varied in the Z direction by removing and adding material to the fracture profile. The creation of a z-varied model geometry is shown in Fig. 5 with a detailed process diagram. The mesh type adopted in this study was predominately quadrilateral based. All models were subject to finer element sizing within the fracture rather than the reservoir zone, making it possible to capture the fluid flow interaction between the two. The z-varied model required a hybrid mesh of quadrilateral and tetrahedral elements due to the nonsymmetrical nature of the geometry (Fig. 6). The total number of elements and nodes in the mesh are 448,000 and 525,000 respectively with the skewness and aspect ratio of 0.1 and 1. Mesh convergence was required to demonstrate the most reasonable element size for the fracture profiles that would ensure an efficient computational time and low discretisation error. The mesh sensitivity study was carried out by simulating the parallel plate fracture profile using 5 Pa pressure inlet condition, 0.2 mD permeability and porous Average Aperture Height=0.35 mm Tortuosity=88%

251
interface condition, shown in Tab. 1, to determine the optimal element size. The results from Tab. 1 suggests that a grid-size of 0.009 produced minimal deviation in maximum pressure and velocity values. Therefore, an element size of 0.009 was selected in this study.  The standard implicit pressure based steady-state solver with double-precision was used as the fluid used for modelling is water (incompressible flow), neglecting any gravitational effects. The laminar viscous model with single phase water was used due to flow resulting in a low Reynolds number. To represent the driving differential pressure forces of reservoir flow more accurately, pressure inlet and pressure outlet boundary conditions for all the models were used as shown in Fig. 6. The momentum equation was solved using the second-order upwind scheme. The interface between the fracture and matrix was varied between the wall condition and the porous condition. When a wall interface condition was used, no fluid can flow from fracture to matrix, and conversely, when the porous condition was used, the fluid can leak off between the fracture and reservoir. A total of 36 simulations were performed for different pressure inlets, permeabilities and fracture boundary conditions. The porosity was assumed to be 20% based on sandstone reservoir. Tab. 2 shows the varying parameters and its range studied in the simulation model.

Results and discussion
The single-phase fluid flow of water through three fracture models namely coarse, parallel plate and z-varied geometry was investigated using the effect of pressure and isotropic permeability to demonstrate the flow behaviour. In addition, the effect of anisotropic permeability on the coarse model with porous fracture-matrix interface condition was established incorporating the effect of inertial resistance as well as viscous resistance. The results obtained are discussed in the following sections.

Effect of pressure
In order to accurately produce output pressure values representative of the actual fractured reservoir, a series of Iso-planes were set up in nine separate locations along with the three models length of 6.6925 mm. Each Iso-plane used the area weighted average pressure to produce either the fracture zone pressure at the selected Iso-plane location or reservoir pressure. The Iso-plane locations were position at 1/10th increments of the total length (6.6925 mm), as shown in Fig. 7. The single-phase fluid flow of water through a macro scale reservoir with one single fracture and isotropic permeability was investigated with varying inlet pressure.

Comparison of porous and wall fracture-matrix interface condition for parallel plate model
In order to fully understand the influence of the fracture-opening, both the fracture and reservoir pressure drops were plotted together at nine iso-plane locations in Fig. 8. Firstly, from Fig. 8 a higher deviation from the reservoir pressure can be seen, primarily in Iso-plane locations 3, 4 and 5. This is due to the predominant area of the reservoir concerning the smaller fracture cross-sectional area, thereby highlighting the necessity to analyse fracture flow independent of the reservoir. An identical trend was observed at 1000 Pa for both the fracture and reservoir pressure, whereby small differences in pressure distribution could only be seen in the pressure contours. The effect of increasing the pressure from 5 Pa to 1000 Pa using a porous condition produced slightly different results. However, the fracture pressure trend remained the same as the wall condition in Fig. 8, the pressure drops slightly for the porous condition, as shown in Fig. 9, at Iso-plane 4.  , vol.16, no.2, pp.245-266, 2020 Figure 8: Parallel plate pressure drop for reservoir and fracture This is because of the porous condition in which fluid can flow between the fracture and the reservoir; therefore, fluid from the surrounding porous reservoir is drawn into the fracture path of low flow resistance. This owes to the slight decrease in pressure shown for the porous condition compared to the wall condition, as the fluid flow is slightly higher due to the reaction between the reservoir and fracture. In Fig. 9 a 3.833% difference between the porous and wall condition pressure values occurs at 5 Pa. To see if the pressure decreased for porous at higher inlet pressures, the same pressure drop was plotted for 1000 Pa. Similarly, the pressure was found to decrease slightly for the porous condition. However, the percentage difference between the wall and porous was much smaller. It can be seen that at Iso-plane 4 the percentage difference between wall and porous for 1000 Pa was 2.336%. It was expected that the percentage difference would be more significant for 1000 Pa due to a higher pressure inducing a higher flow from the porous matrix to the fracture opening, thus causing a larger drop in pressure. In comparing Figs. 9 and 10, it is clear that this only occurs in Iso-planes 6-8, or further along the fracture where the flow has become established in the reservoir, and the effect of the higher pressure on the flow from the porous matrix to the fracture can be felt.

Comparison among parallel plate, coarse and z-varied models
To investigate the effect of fracture coarseness on the coarse and z-varied model, it was imperative to compare the pressure plots and resulting in maximum velocity to the parallel plate model. According to Rasouli et al. [Rasouli and Rasouli (2012)], one of the main differences with a coarse fracture in comparison to a parallel plate is an increase in permeability and therefore, flow rate due to rougher surfaces exhibiting larger mean apertures. To examine this theory, the fracture pressure drop was plotted for both the coarse and parallel plate model, as shown in Fig. 11. The pressure at iso-plane 5 shows that the parallel plate pressure dominates at 3.96 Pa compared to 3.23 Pa for the coarse model. Further down the fracture, at 3.6 mm the pressure relationship changes, the coarse fracture pressure now dominates, particularly at iso-plane 8 where the coarse pressure is 2.18 Pa compared to 1.28 Pa for the parallel plate. The location of the minimum aperture for the coarse and parallel plate models differ, as shown in Fig. 12. Therefore, it was expected that the parallel plate maximum velocity would be higher due to the lower minimum aperture of 0.0316 mm compared to 0.04987 mm for the coarse model. However, this contradicts the theory identified by Rasouli et al. [Rasouli and Rasouli (2012)] that the coarse model should give rise to a higher velocity compared with the parallel. According to Fig. 12, the location of the minimum aperture for the parallel plate model lies between Iso-planes 5-6 while the minimum aperture for the coarse model lies between Iso-plane 8-9. Taking the gradient of these zones shows that the parallel model has a steeper pressure drop, thus proving that the parallel plate model should in theory exhibit a higher maximum velocity than the coarse model. Through fluid modelling it was discovered that the coarse model generated a higher maximum velocity than the parallel plate model, thus agreeing with Rasouli et al. [Rasouli and Rasouli (2012)]. Tab. 3 shows the maximum velocity for the coarse model, which was found to remain high in relation to the parallel plate model for both 5 Pa and 1000 Pa at both porous and wall interface condition. Furthermore, on comparison of the z-varied model with the parallel plate model and coarse model in Fig. 13 shows that between 0.001 m and 0.004 m the z-varied model exhibits a higher pressure compared to that of the parallel plate model. After 0.004 m the z-varied model exhibits a lower pressure compared with the parallel plate. The pressure dropping to lower than the parallel plate near the end of the fracture is due to the z-varied model possessing a lower surface area due to the amount of material removed. The flow in the fracture will be induced due to the removal of the material causing the velocity to increase and therefore, the pressure to decrease. However, this does not justify why the pressure appears higher in the early stages of the fracture for the z-varied model in comparison to the parallel plate, shown in Fig. 13. Therefore, it can be said that the effect of the reduced material is only felt towards the end of the fracture where the fracture fluid has built up speed from the reduced area earlier in the fracture. It can be said that the flow is developing, and slowly responding to the effect of the lowered fracture opening volume.

Figure 13:
Pressure drop for coarse, z-varied and parallel at 0.2 mD 5 Pa At a higher inlet pressure of 1000 Pa the pressure drop was identical to that of Fig. 13, showing that a higher pressure has little to no effect on the pressure distribution of the zvaried model. The investigation covers the effect of pressure on three models; parallel plate model, coarse model and a parallel plate model with apertures varied in the z-direction. Close attention to the distribution of pressure and velocity of the flow for all the models were made with the comparison. Each model has been investigated for the effect of adding a porous condition across the fracture reservoir interface, allowing both the reservoir flow and fracture flow to interact with varying inlet pressures. Comparing all the three models, it can be seen from Tab. 3 that the z-varied model yields a higher maximum velocity at 1000 Pa for the porous condition. At a low permeability of 0.2 mD and a high inlet pressure of 1000 Pa for the porous condition, more fluid is drawn into the matrix from the porous matrix. Therefore, a higher maximum velocity occurs. The z-varied model, due to the reduced area means that the already increased velocity at this condition is further increased due to the removal of material. This shows why the z-varied model maximum velocity appears higher in relation to the equivalent parallel plate and coarse models at 1000 Pa, 0.2 mD for the porous condition. At the wall condition for 1000 Pa, 0.2 mD the z-varied model appears to have the lower maximum velocity in comparison to the coarse and parallel plate model. This shows the importance of a permeable wall in influencing the flow rate in the fracture zone. Not only is the permeability of the fracture important but the pressure needs to be sufficient enough to influence the rate at which the flow is drawn into the fracture. This is why the maximum velocity is smaller for the zvaried model at a low inlet pressure of 5 Pa for both the porous and wall conditions in comparison to the parallel and coarse models.
In order to verify that the porous condition yields a slightly lower pressure in comparison to the wall condition, the pressure drop was plotted for all the models with wall and porous condition. It is evident from Fig. 14 that the trend is the same in terms of the porous model producing a slightly lower pressure in comparison to the equivalent wall condition model. It is also clear that the pressure drop is more significant between the porous and wall condition for the coarse model, indicating that the coarser geometry yields the higher maximum velocity for the 5 Pa, 0.2 mD model. This can be confirmed by Tab. 3, where the coarse porous model for 5 Pa, 0.2 mD produces the highest maximum velocity of 0.00318 m/s out of the corresponding coarse and parallel plate models.
The results of the present study were compared against the research study by Crandall et al. [Crandall, Ahmadi and Smith (2010)]. The simulation was performed with a similar geometry, mesh, input boundary conditions and modelling parameters. The steady-state laminar viscous model with single phase water was used to simulate the fluid flow. The simulation parameters used are as explained in Section 2.2.

Figure 14:
Pressure drop for all the models at 5 Pa, 0.2 mD (wall and porous condition) Fig. 15 shows the pressure contour plot of the current study for the fracture and matrix, and Fig. 16 shows the comparison of pressure drop profile in the fracture along the nondimensional fracture length, for the inlet pressure, =1000 Pa and permeability=2000 mD.
The results show a good match with the research study [Crandall, Ahmadi and Smith (2010)]. A slight deviation in the results can be attributed to the difference in the fracture geometries concerning aperture. As mentioned earlier, the geometry used in the current study is originally from an experiment from Karpyn et al. [Karpyn, Alakmi, Parada et al. (2003)], and the same geometry is used by Crandall et al. [Crandall, Ahmadi and Smith (2010)] in their study. In the current study, the image was digitised using plot digitiser software for the analysis, which can result in a slight difference in the geometries.

Effect of permeability 3.2.1 Parallel plate model
Permeability is the rock property and describes the ability for fluids to flow through rocks. It is also an indication of connectivity between pores in a reservoir, the higher the connectivity, the lower the resistance to the transfer of fluids in the reservoir. In ANSYS Fluent pathlines can be used to analyse the flow of particles through the pores of a reservoir. Tab. 4 demonstrates the different particle distributions generated for the various pressures and permeability's, under porous interface condition. Particle data (number escaped and number trapped) was measured for 1000 steps at a path skip of 200 and shown in Tab. 4. According to Tab. 4, the porous interface condition for parallel plate indicates that fluid is being interacted between the reservoir and fracture. At 5 Pa when the permeability was increased from 0.2 mD to 2000 mD the particle tracking highlighted that more flow could transport through the porous media, due to improved connectivity, as expected. For 5 Pa, 0.2 mD there appears to be a higher number of particles of 67.9% being trapped compared with 5 Pa, 2000 mD, where 21.1% of particles become trapped. Due to the interaction associated with the porous interface condition, many fluid particles are pursuing the path of least resistance and in doing so are becoming trapped in numerous eddies formed in the porous matrix. This difference can be due to the ease of flow through the porous matrix, higher pore connectivity from the higher permeability induces more flow. Tab. 4 suggests that when the fracture wall is porous and permeable a higher inlet pressure of 1000 Pa at 0.2 mD dramatically improves the flow through the porous media and the path line plot shows enhanced flow in the porous media. According to Tab. 4, 43.9% of particles were trapped for the 1000 Pa, 0.2 mD condition, compared with a higher 67.9% particles being trapped in the 5 Pa, 0.2 mD condition. This indicates that for a porous interface condition the higher-pressure inlet of 1000 Pa has a much more significant effect on the rate at which particles are trapped, compared with wall condition that theoretically has no fluid particles trapped because of mass conservation.  In order to show the effect of increasing permeability on the pressure drop for the fracture and reservoir, the pressure drop was plotted for 1000 Pa (Fig. 17). According to Tab. 5, at 2000 mD the pressure can be seen to drop slightly compared with that of the 0.2 mD model, this can be seen to occur early on in both the fracture and reservoir zone. At a low permeability, higher pressure is maintained more effectively than at a greater permeability. However, this trend switches round at Iso-plane 6 where the lower permeability is shown to have a higher pressure. This can be attributed to the minimum aperture, which has drastically increased the velocity and correspondingly reduced the pressure. At a lower permeability, the pressure takes longer to recover, whereas the high permeability allows for the flow path to improve and recover quicker. The lower permeability of 0.2 mD means there is a limited area for the pressure to recover in the area increase following the minimum aperture. At 2000 mD there is a higher porous area for the pressure to recover in.

Comparison among parallel plate, coarse and z-varied models
Tab. 6 shows the comparison of the parallel plate, coarse and z-varied models, based on particle trapped data from the path line plots. The data from Tab. 6 was conveyed via the bar chart presented in Fig. 18. It is clear that the z-varied model yields the highest percentage of trapped particles of 69.23% at 5 Pa and 0.2 mD. This can be attributed to the removal of material resulting in more zones of higher velocity flow from an apparent reduction in fracture area. More zones of narrow fracture aperture lead to a higher formation of eddies due to more flow encountering fracture expansion from the narrow apertures. Zones of very low pressure occur in a region of slightly higher-pressure due to the rapid change in the area which has caused the velocity to increase and thus the pressure to decrease. The eddies that form in these low-pressure zones draw flow particles in and cause them to become trapped, hence the higher number of trapped particles in the z-varied model. According to Fig. 18, this phenomenon only occurs at a low pressure of 5 Pa, when the fluid has insufficient momentum combined with a low fluid conductivity to overcome the forces drawing the flow particles into the eddies.

Effect of anisotropic permeability-coarse model
This paper has explored the effect of varying isotropic permeability, however, in a realistic reservoir, the permeability will differ in the x, y and z planes. The analysis is done for the anisotropic condition in the coarse porous model. From Tab. 7 the three permeability configurations, A21, A23 and A32 were analysed for 5 Pa and 1000 Pa and effectively compared with the isotropic A22 model. When the permeability in the x-direction was set to 2000 mD, the main difference encountered was the lower pressure distribution towards the end of the fractured reservoir Section. The difference occurred between 0.005 m to 0.006 m where the pressure for Kx=2000 mD was 10.4% lower than the pressure for Kx=0.2 mD. This can be confirmed by Fig. 19 where the pressure drop trend for the three permeability conditions is almost identical, and the higher x-direction permeability yields the lower pressure due to a higher overall flow rate in the fractured reservoir. The effect of the higher x-direction permeability can be seen to only affect the fracture pressure at around 0.003 m. At this point, the effect of the higher formation permeability in the x-direction can be felt and continues to affect the fracture flow as the length of the fractured reservoir increases. This suggests that the fractured reservoir sample in question is too small to show the effect of a higher formation permeability in the x-direction. 2 3 x-y 1 2 A21 A23 3 A32 A21 => KX=0.2 mD, Ky=0.2 mD, Kz=1000 mD A23 => KX=2000 mD, Ky=2000 mD, Kz=1000 mD A32 => KX=1000 mD, Ky=1000 mD, Kz=2000 mD  Towards the end of the fractured reservoir model, the effect of the higher x-direction permeability can be seen to affect the distribution of the lower pressure, as shown by the pressure contour in Fig. 20. For the A23 permeability arrangement, the pressure was probed after Iso-plane 9 and showed that a percentage difference of 32.1% existed compared to that of the isotropic A22 model. Therefore, it was established that an even greater percentage difference exists when comparing A23 to the isotropic A22 model. The smaller 10.4% difference between the pressure for the A23 and A21 permeability arrangements confirms that the higher permeability of 1000 mD in the z-direction influences the area in which fluid can flow and thus improves flow through the matrix reducing the corresponding pressure difference.

Conclusions
The main findings from the numerical modelling of single-phase fluid flow through complex fractured reservoirs are as follows: • The fluid flow through fractures is highly dependent on the fracture aperture.
• The nature of the fracture-matrix interface condition can play a significant role in fluid flow through fractures. At the same pressure and permeability conditions, the porous condition was found to reduce the pressure by 2-10% in comparison to the wall condition. This shows the importance of a permeable wall in influencing the flow rate in the fracture zone. Not only is the permeability of the fracture important but the pressure needs to be sufficient enough to influence the rate at which the flow is drawn into the fracture. • The z-varied model provided the more representative with the parallel plate, and coarse geometries showed the highest number of particles being tracked at low pressure, low permeability. This can be attributed to the higher formation of eddies caused by a higher number of narrow apertures in the z-varied model.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.