Fluid flow through porous media using distinct element based numerical method

Many analytical and numerical methods have been developed to describe and analyse fluid flow through the reservoir’s porous media. The medium considered by most of these models is continuum based homogeneous media. But if the formation is not homogenous or if there is some discontinuity in the formation, most of these models become very complex and their solutions lose their accuracy, especially when the shape or reservoir geometry and boundary conditions are complex. In this paper, distinct element method (DEM) is used to simulate fluid flow in porous media. The DEM method is independent of the initial and boundary conditions, as well as reservoir geometry and discontinuity. The DEM based model proposed in this study is appeared to be unique in nature with capability to be used for any reservoir with higher degrees of complexity associated with the shape and geometry of its porous media, conditions of fluid flow, as well as initial and boundary conditions. This model has first been developed by Itasca Consulting Company and is further improved in this paper. Since the release of the model by Itasca, it has not been validated for fluid flow application in porous media, especially in case of petroleum reservoir. In this paper, two scenarios of linear and radial fluid flow in a finite reservoir are considered. Analytical models for these two cases are developed to set a benchmark for the comparison of simulation data. It is demonstrated that the simulation results are in good agreement with analytical results. Another major improvement in the model is using the servo controlled walls instead of particles to introduce tectonic stresses on the formation to simulate more realistic situations. The proposed model is then used to analyse fluid flow and pressure behaviour for hydraulically induced fractured and naturally fractured reservoir to justify the potential application of the model.

Abstract Many analytical and numerical methods have been developed to describe and analyse fluid flow through the reservoir's porous media. The medium considered by most of these models is continuum based homogeneous media. But if the formation is not homogenous or if there is some discontinuity in the formation, most of these models become very complex and their solutions lose their accuracy, especially when the shape or reservoir geometry and boundary conditions are complex. In this paper, distinct element method (DEM) is used to simulate fluid flow in porous media. The DEM method is independent of the initial and boundary conditions, as well as reservoir geometry and discontinuity. The DEM based model proposed in this study is appeared to be unique in nature with capability to be used for any reservoir with higher degrees of complexity associated with the shape and geometry of its porous media, conditions of fluid flow, as well as initial and boundary conditions. This model has first been developed by Itasca Consulting Company and is further improved in this paper. Since the release of the model by Itasca, it has not been validated for fluid flow application in porous media, especially in case of petroleum reservoir. In this paper, two scenarios of linear and radial fluid flow in a finite reservoir are considered. Analytical models for these two cases are developed to set a benchmark for the comparison of simulation data. It is demonstrated that the simulation results are in good agreement with analytical results. Another major improvement in the model is using the servo controlled walls instead of particles to introduce tectonic stresses on the formation to simulate more realistic situations. The proposed model is then used to analyse fluid flow and pressure behaviour for hydraulically induced fractured and naturally fractured reservoir to justify the potential application of the model.

Introduction
Fluid flow through porous media has been the subject of interest in many areas, such as petroleum and resource engineering, geothermal energy extraction, and/or ground water hydrology, etc., for many years. Many analytical as well as numerical models have been developed to explain fluid flow through porous media. Although most of these models work well with reasonable accuracy in the case of reservoir consisting of homogenous media with simple geometry, they encourage huge uncertainty with erroneous results in the case of heterogeneous and discontinuous media, especially in presence of natural fractures and interacted induced hydraulic fractures with complex geometry. In addition, numerical or analytical methods of solutions depend on the geometry of the reservoir, fluid flow type (Linear, Radial and Spherical), fluid flow regime (Transient, Late Transient, Steady State, Semi-Steady state), as well as discontinuity (conductive discontinuity with different permeability, sealed discontinuity with low permeability). To combine all of these factors to get a solution that can describe the flow, many assumptions need to be made. The more assumptions that are integrated into the solution, the more susceptible the solution will be to incur errors in results. Another problem with these methods is that they are not unique solutions. At any time, if one of the factors that are mentioned earlier is changed, the whole solution might change. A robust model is always desirable to address all of these issues and provide better and more accurate solution. In this view, study has been focused to develop a numerical tool to analyse fluid flow in the above mentioned complex scenarios.
Particle flow code (PFC) developed by Itasca consulting group which is based on distinct element method (DEM), is considered as a numerical tool to simulate fluid flow in porous media in this study. The initial fluid flow model was developed by Itasca. Further modifications were made on the model to simulate more realistic situations and validate the simulation results. The method developed is independent of the reservoir geometry, discontinuity, fluid flow type and regime; and is found to be more appropriate to simulate the reservoir that is heterogeneous in nature in relation to both the media and complex geometry (grain shape, size, as well as pore geometry).
A numerical model to simulate the fluid flow through heterogeneous porous media, especially in presence of natural fractures interacted with hydraulic fracture is presented in this paper. In this study, two cases of laminar and radial fluid flow conditions are simulated using the numerical model developed based on distinct element method. The accuracy of the model is validated by comparing the simulation results with analytical results.
After validation, this model is used to simulate fluid flow in a reservoir that is hydraulically fractured and contains two sets of natural fractures. Based on the results, it is demonstrated that proposed DEM based model can potentially be used to analyse fluid flow through complex reservoirs.

Discrete element method
Reservoir rock consists of grains, pores which are filled with pore fluids and possibly joints and faults or in a general term discontinuity which may or may not be filled with cement. If the size of discontinuities is not in the scale of the reservoir rock, a continuum based model may be accurate enough for simulation by incorporating some modification in the model to include the effect of these features. However, if the size of discontinuities are comparable with the size of the rock, continuum based models may lose their accuracy. In this case, a discontinues based model will provide better results (Morris et al. 2003). One of these discontinues models is Discrete Element Method, which is a family of numerical methods that defines the domain as a combination of independent elements. This method is mostly used for granular media, fractured rock systems, systems composed of multiple bodies in mechanical engineering and also fluid mechanics. During 1970s-1980s, this method has rapidly developed for geological and engineering applications. The major breakthrough was by Cundall in 1971, to study rock mechanical properties. In 1979, Cundall and Stark applied this method for soil mechanics. In Discrete Element Method, the independent particles can be either rigid or deformable, and can have circular or polygonal shape (Jing and Stephansson 2007). A great advantage of Discrete Element Method compared to continuum based model is that meshes are built by individual elements and there is no need for remeshing as the simulation progresses. Figure 1   Distinct element programs have been developed based on Distinct Element Method (DEM) which is a sub-classification of Discrete Element Method. In this method, contacts are deformable and discretised elements can be either rigid or deformable (Cundall and Hart 1992). The solution scheme is based on explicit time stepping which time steps are chosen so small that the disturbances introduced by a single particle cannot propagate beyond neighbouring particles (Cundall and Strack 1979). Detailed description of the method can be found in (Cundall 1988) and (Hart et al. 1988).
PFC2d (Particle Flow Code in two dimensions) is a DEM based commercial software developed by Itasca Consulting Group. In this software, discretised bodies are composed of rigid circular particles that can have a random distribution of radius size from a range defined by the user. The analysis is based on Force-Displacement calculation for individual particles and applying Newton's second law for calculating velocity and position of particles in each time step (Itasca 2008a). Figure 3 depicts the general algorithm used in PFC: Figure 4 shows a collection of particles that have been generated in PFC. Each particle in contact with other particles can cause normal and tangential forces. The magnitude of these forces depends on the overlap of particles, elastic properties of particles, contact model and contact model properties.
More information about the formulation and analysis procedure can be found in PFC2D manual (Itasca 2008a).

DEM fluid flow
In PFC2D, porous medium, through which fluid flows, is considered to be composed of individual circular particles, which are connected together by contact or parallel bond. The void space between particles is assumed to be filled with fluid, which flows between these void spaces. To characterize fluid flow between these void spaces, it is required to define the domain term. A domain is defined as a closed loop polygon by particles that are connected to    (Huynh 2014) show boundaries of domains and they connect centre of particles that build the domain.
Fluid flow happens between domains through a pipe centred at the contact point between each two particles. Pipe length is the sum of two particles radius. Aperture between parallel plates is denoted as ''w''. The depth of the pipe is equal to unity (Itasca 2008b). Figure 7a shows two domains connected by the pipe. Domains are denoted as Domain i and Domain j , and the pipe connecting them is shown as a rectangle in red colour. Figure 7b shows the pipe. Aperture of the pipe is w and its length is L P . Depth of the pipe which is in out of plane direction is equal to one.
Fluid flow through pipe is governed by the Poiseuille Equation for fluid flow through parallel plates (Garde 1997): where; q: Flow Rate, P i : Pressure in Domain i, P j : Pressure in Domain j, l: Fluid Viscosity, L p : Pipe Length.
Pipes can be defined between two particles, only if they are in contact. However, after they have been initialized they will exist even if particles detach from each other. When particles are in contact, the aperture of the pipe will be zero. But to take into account the macroscopic permeability of the rock and to overcome the 2D limitation of the simulation, ''w'' will be set to a number greater than zero (Itasca 2008b).
After setting the initial value for aperture, its value should also take into account the nature of the contact between its two adjacent particles. That means, if particles are still in contact and apply a compressive force on each other, ''w'' can be found by Eq. 2 (Itasca 2008b): w 0 is the initial aperture. F 0 is fixed value and is the amount of normal force that changes the aperture to half of its initial value. F is the compressive force between particles and its value can change. If the value of F is much lower than F 0 , the value of w would not change appreciably. If the value of F is negative (i.e. particles bond is under tension), the aperture value is obtained by Eq. 3 (Itasca 2008b): In this equation, g denotes the gap or the distance between particles. m is a calibration constant and can have a value between zero to one. Macroscopic permeability value of the rock can be reproduced by adjusting values of w 0 , F 0 and m. Each domain has pressure communication with other surrounding domains. Total flow in and out of each domain is P n i¼1 q i with ''n'' being the number of surrounding domains that is different for each individual domain with a minimum value of one and maximum value can be any number greater than zero. Total flow volume into a domain in one time step is given by Eq. 4 (Itasca 2008b): ''n'' is different for different domains and is equal to number of surrounding domains and Dt is time step. In each time step, domains experience a mechanical volume change because of movement of particles. These movements are because of changes that occur in forces between particles. Total pressure change in one time step for a single domain is (Itasca 2008b): DV d is mechanical volume change of the domain, V d is the volume of the domain and K f is the bulk modulus of the fluid. After a sample is generated, it will be enclosed by four frictionless plates as shown in Fig. 8. These plates can move independent of each other toward or away from sample. Plates have no interaction on each other and only interact with particles of the sample. The purpose of these plates is to introduce principle stresses on the specimen to resemble tectonic stresses that are present underground.
Because of principle stresses, particles exert normal and tangential forces on each other. In addition to these forces, fluid inside pore spaces also exerts some force on particles. Figure 9a shows a domain with pressure P i and Fig. 9b shows particle 1 and forces on it that are generated because of pressure in domain. Resultant of these forces is F.
The magnitude of force that is exerted on Particle 1 by fluid pressure in domain i is the product of pressure by the area of the particle that is exposed to domain i.
The depth of the particle is equal to 1. So the area is equal to 1 multiplied by the length of the arc. Length of arc is equal to radius of the particle multiplied by the angle that forms between the lines joining the centre of the particle to its neighbouring particles. Direction of force is outward from the domain in the direction of the line that divides the arc into two halves.
If h in Eq. 6 is greater than p, then it should be subtracted from 2p and the result be substituted instead of h. This is shown in Fig. 10. It is evident from figure that forces that are applied to the shaded section of the particle will balance each other out. The net force that remains on the particle as a result of pressure in domain i will be equal to the forces that are applied on the arc between the dashed lines. Also from figure, it can be seen that For simplicity of calculations, the mechanical volume change in domains is neglected as its value is very small and will not make a noticeable change in results. On the other hand, domain volumes will be updated in every time step. The solution to fluid flow alternates between flow through pipes and pressure adjustments between domains. This means, in each time step, the fluid flow through pipes is calculated and the total net flow to or from each domain will cause pressure change in domains. For stability Fig. 7 a Pipe connecting two domains. Pipe is shown in red colour. b Pipe with length L P , width w, and depth 1 Fig. 8 Sample with principle stresses acting on its sides analysis, a critical time step needs to be calculated. The procedure is to first calculate the critical time step for each domain and then take the minimum time step of all critical time steps as the global time step (Itasca 2008b).
In every time step, total flow into a domain because of pressure perturbation, DP p is calculated by Eq. 8 (Zhao 2010).
N is the number of pipes for each domain, and R is the average of radius of particles that surround the domain. Solving for DP p : The pressure disturbance applied to domain because of the above flow rate is (Zhao 2010): To confirm stability, the pressure response needs to be less than or equal to pressure perturbation: If we replace Eqs. 9 and 10 in Eq. 11 we will get Eq. 12: Solving Eq. 12 for Dt will result in Eq. 13: To simplify the right hand side of Eq. 13, DV d Q is neglected and the remaining terms are multiplied by a safety factor (Zhao 2010): The fluid flow calculation is explicit in time. Figure  It can be concluded from this section that the fluid flow at microscopic level is independent of the reservoir geometry because the fluid flows between the domains, can be created for any reservoir shape. This characteristic makes this method applicable for any reservoir geometry. Also, this method is independent of flow type (i.e., linear, radial, etc.) and flow regime (i.e., transient, late transient, steady state or semi steady state), and can be applied for any flow type and regime.
A discontinuity such as a fracture or a joint has a permeability that can be different from matrix permeability. These discontinuities can be incorporated into the system Fig. 9 a Domain with pressure P, b Pressure applied to part of particle 1 that is exposed to domain 1 and generated force F because of pressure P by using smooth-joint model. Domain pipes that are on smooth-joint can have different aperture (w) and by adjusting these apertures, discontinuity permeability can be reproduced. In this way, the only difference between pipes that represent discontinuity and pipes that represent pore throat is their aperture. Following this method simplifies the incorporation of discontinuities without requiring developing a whole new system to describe fluid flow. Analytical models, on the other hand, may require development of whole new solution as a new type of discontinuity is presented in the system. Figure 12 illustrates a sample with two sets of joints with dip angles of (60°) and (-20°).

Analytical methods
To validate numerical model, two cases of linear fluid flow and radial fluid flow in porous medium is considered. The formation in both cases is finite. For each case, the analytical formula with its boundary and initial conditions is presented in Sects. ''Linear fluid flow condition'' and ''Radial fluid flow condition''. Derivation of these analytical equations is presented in Appendices A and B. Solutions of both cases are used in the comparison section to compare the results of numerical model against analytical models.

Linear fluid flow condition
One dimensional fluid flow is considered in a sample with dimensions of L Â H Â W. The initial pore pressure of the sample is set equal to P i . Boundaries have constant pressure. The pressure at one end of the sample is P 1 and on the other end is P 2 . Equation 15 shows the initial condition. Equations 16 and 17 show the boundary conditions: Equation 18 shows the linear form of pressure diffusion equation in porous media (Donnez 2012): a 2 is called hydraulic diffusivity and is equal to fluid mobility divided by fluid storability (Donnez 2012): The solution of Eq. 18 is presented in Appendix A. The final solution is shown in Eq. 20: To simplify calculations, the concept of dimensionless parameters is used. Dimensionless pressure, position and time are defined as (Ahmed and McKinney 2011): If P 1 is set equal to P i and P 2 is set equal to zero, after re-arranging Eq. 20 and using Eqs. 21-23 dimensionless pressure is: Equation 24 is the dimensionless form of pressure diffusion equation of laminar fluid flow in a sample with initial pore pressure of P i and constant boundary pressures of P 1 ¼ P i and P 2 ¼ 0.

Radial fluid flow condition
Radial fluid flow is considered for a sample with external radius (R e ) and wellbore radius ðR w Þ. Initial pore pressure is P i . Pressure at outer boundary ðR e Þ is kept constant at P i while wellbore pressure changes to ensure constant wellbore flow rate. Equations 25 and 26 show inner and outer boundary conditions (Ahmed and McKinney 2011). Equation 27 is the initial condition.
The radial form of pressure diffusion equation in porous media is shown in Eq. 28 (Dake 1978): Solution of Eq. 28 is shown in Appendix B. The final solution in dimensionless form is: Equation 29 is the dimensionless form of pressure distribution across a circular reservoir with constant outer boundary pressure and constant wellbore flow rate. k n in this equation are roots of Eq. 30. Y 1 is first order of second kind Bessel function, Y 0 is zeroth order of second kind Bessel function, J 0 is zeroth order of first kind Bessel function and J 1 is first order of first kind Bessel function.
Knowing R De , Eq. 30 can be solved to get the values of k n . There are infinite numbers of k n that will satisfy this equation and they are called eigenvalues of this equation. Dimensionless parameters in Eq. 29 are defined as: R D is dimensionless radius, t D is dimensionless time and P D is dimensionless pressure.

Comparison of numerical and analytical models
Linear fluid flow with constant boundary pressures Figure 13 shows a sample before setting the pore pressure (a) and the same sample with pore pressure being set (b). Right hand side of the sample in (b) has no pore pressure as its pore pressure is set equal to zero and will be kept at zero during fluid flow. The left hand side pressure will be kept constant at initial pressure. Sample length is 14 meter, its height is 14 meter and its width which is an out of plane dimension is equal to one meter. The sample length between the left boundary and the right boundary is 13.5 meters. Brown circles in Fig. 13b represent domain pressure. Figure 14 shows the simulation results at four different times. Each point shows the pressure of a domain. Vertical axis shows pressure in Mega Pascal (MPa) and horizontal axes show the x and y position of the domain. The pressure at left hand side is kept constant at initial pressure and pressure at right hand side is kept constant at zero. At the beginning of simulation, flow rate on the left hand side is zero as pressure reduction wave has not reached it and the flow rate on the right hand side is 1:101 Â 10 À3 m 3 =s. The flow rate quickly drops on the right hand side as the pressure on the right hand side start to fall down. As soon as the pressure reduction wave reaches the left hand side, the flow rate on left hand side start to increase. Flow rate keeps increasing on left and falling on right until steadystate flow rate is established. At steady-state condition, flow rate on both sides is equal to 1:365 Â 10 À5 m 3 =s. Figure 15 shows pressure distribution across sample with each point representing the pressure of its domain. At t ¼ 15067 s, a steady state flow regime is established.
To make sure that simulation results are correct, they are compared against analytical results. To do so, data in Fig. 15 are converted to dimensionless form by using Eqs. 21-23.
The permeability of the simulated sample can be obtained by using Darcy equation for linear flow in steadystate condition as given by Eq. 34.
By inserting the values of different parameters from Table 1 into Eq. 34, permeability is found to be 2:67 Â 10 3 md. To confirm that this is the correct value, another simulation has been conducted with different initial pressure of 8:00 MPa. Steady state flow rate was 2:18 Â 10 À5 m 3 =s. Rest of the parameters were kept constant. By inserting new values for initial pressure and flow rate, permeability is calculated to be 2:67 Â 10 3 md, which is same as the value calculated before. Dimensionless times and position are also inserted in Eq. 24 to compare simulated versus analytical results. Equation 35 is the dimensionless time in field units. Table 2 shows simulation time in seconds and day and calculated dimensionless time.
The units for different parameters are: k: md, t: day, l: cp, c t : psi -1 , L: ft. Figure 16 shows results of simulation in dimensionless form versus analytical results. It shows that data from simulation match very well with analytical results and validates the model for linear fluid flow.

Radial fluid flow
Radial fluid flow is simulated for a reservoir with constant external boundary pressure and constant wellbore flow rate. Figure 17a shows the reservoir with wellbore in the centre.  Figure 17b shows the reservoir with boundary pressure set. Pressure at distances more than external radius R e is constant. Wellbore radius is R w and height of the reservoir which is out of plane dimension is equal to one. Each brown circle represents pore pressure of a domain which is kept constant at 7 MPa at distances greater than external radius. Figure 18 shows the reservoir with initial pore pressure set to 7 MPa. The initial and boundary pressure are equal. Figure 19 shows simulation results within the external radius of the reservoir. Before starting production from wellbore, the pressure across the whole reservoir is same and equal to P i . As production starts with constant wellbore flow rate, pressure starts to decrease near the wellbore. As simulation continuous, the pressure reduction wave propagates toward the outer boundary. Before pressure reduction arrives at outer boundary, flow rate at outer boundary is zero. As soon as the pressure reduces near the outer boundary, the flow rate starts to increase. Flow rate keeps rising until steady-state condition is reached. At steadystate condition, outer boundary flow rate is equal to wellbore flow rate. Table 3 shows simulation parameters: At steady state condition, wellbore pressure becomes constant and equal to 3:49 Â 10 6 Pa and is used to determine the permeability of the reservoir. Steady state radial flow in field units is: The units for different parameters are: k: md, P: psia, l: cp, h: ft, R: ft. Using Eq. 36, permeability is found to be equal to 2:44 Â 10 3 md. Figure 20 shows simulation results of pressure versus radius. Vertical axis shows pressure in Pa and horizontal axis shows radius in meters. Each point shows the pressure of its domain versus the radial distance of the domain with respect to wellbore centre. The figure shows that pressure at outer boundary is kept constant while wellbore pressure keeps reducing until steady state is reached. To make sure that simulated results are correct, they are compared against analytical results. To do so, Eq. 29 is used to find the dimensionless pressure versus dimensionless radius. In this equation, dimensionless time is required which is calculated by Eq. 37: The units for different parameters in Eq. 37 are: k: md, t: hr, l: cp, c t : psi -1 , r w : ft. Table 4 shows simulation time in the left column and calculated dimensionless time in the right column.
Equation 29 also requires values of k n , which are the roots of Eq. 30. A plot of ½Y 1 k ð ÞJ 0 kR De ð ÞÀJ 1 k ð ÞY 0 kR De ð Þ versus k is drawn in Maple and is shown in Fig. 21. R De is equal to 14.26. It is evident from the graph that values of function approaches zero very quickly as value of k increases. So to solve Eq. 29, only first few roots of the function will be sufficient to get acceptable results. This equation is solved in maple for first 50 roots. The values of k are shown in Table 5: Simulation and analytical dimensionless pressures versus dimensionless radius are calculated and plotted in Fig. 22. This plot shows that simulation results match with analytical results very well and validates the applicability of the model for radial fluid flow.

Fractured reservoir
In this section, fluid flow in a fractured reservoir is simulated. Figure 23 shows the reservoir with wellbore at the centre, a bi-wing hydraulic fracture and two sets of natural fractures. The angle of the hydraulic fracture is 0 and first and second sets of natural fractures have 65 and À45 angle with respect to x axis, respectively. Each brown circle shows the domain pressure at distances greater than  external radius of the reservoir. Pore pressure of the reservoir is not shown in this picture, but its initial value is same as the boundary pressure and equal to 7 MPa. Boundary pressure is kept constant and reservoir pore pressure is allowed to drop to keep a constant wellbore flow rate at 2 Â 10 À10 m 3 =s. Permeability of rock matrix is 0:458 md and fracture permeability is 2:67 Â 10 3 md. Figure 24 shows two-dimensional view of pressure distribution in the reservoir at six different times. In this figure, it can be seen that pressure in the fractures drops at a faster rate with respect to rock matrix pressure. Figure 25 shows a three-dimensional view of the pressure distribution at time 38287:82 s. Because all the fractures have a Tshaped connection at intersection with other fractures, their pressure shows a V-shaped plot and the pressure data of corresponding fracture can be identified easily. Figure 26 shows three-dimensional view of the pressure distribution at six different times. In Fig. 26a and b, only pressure drop in hydraulic fracture can be seen. As the pressure reduction wave arrives at the intersection of hydraulic fracture and first set of natural fracture, pressure start to drop in natural fracture and this can be seen in Fig. 26c. Figure 26d shows the start of pressure drop in second set natural fractures. The permeability of fractures is adjusted by adjusting the aperture of the pipes that fall on the position of the fracture. For simplicity of observing the clear plot of pressure data in fractures, all fractures permeability was set to be same. But to account for the permeability difference between different sets of natural

Potential application of the model
As stated earlier in introduction, fluid flow through porous media has been the subject of interest in many fields of study. The proposed model is expected to be considered in areas associated with fluid flow through complex heterogeneous porous media subjected to anisotropic stress system. The model is also expected to be used to study many issues related to fluid production from, and injection into porous media with complex situations, especially in case of naturally fractured media, and hydraulic fractures inter-  1. Oil and gas flow in oil and gas reservoirs. 2. Injection of water or surfactants into oil and gas reservoirs. 3. Both the injection and production simulation in oil and gas reservoirs. 4. Single stage as well as multi stage hydraulic fracturing of the oil and gas reservoirs such as shale gas reservoirs. 5. Water flow in mines both in rock matrix as well as in joints, faults and fractures for the application in mining industry. 6. Water flow underground in soil to be used by civil engineers for simulating water flow into tunnels. 7. Water movement in soil to be used by agricultural engineers to simulate the rate of hydration, dehydration or draining of soil.

Conclusion
Deriving analytical expressions to describe fluid flow in porous medium is a complex task. This is because for any change in the reservoir geometry or any change in the condition of fluid flow (e.g., transient, late transient, semi-steady state or steady state) a new analytic expression needs to be developed. In this view, a DEM based numerical model is proposed to analyse the fluid flow through reservoir's porous media with complex characteristics, especially in the case of existence of natural fractures, hydraulic fractures and interaction of hydraulic and natural fractures for any condition of fluid flow. Proposed model is used to simulate and analyse some of these field representative cases as an example case studies. Both simple and complex cases are considered in this study. The simple case is used to validate the accuracy of the model. The DEM model that was used in this study is observed to be independent of reservoir geometry as well as the condition of fluid flow since it was shown that it worked for both linear and radial flow without modification. The simulation results are found to be in good agreement with analytical results. It is also demonstrated that the model can potentially be used as a powerful numerical simulation tool to handle both simple and complex reservoir conditions such as complex formations with irregular shapes, and complex set of discontinuity and fluid flow regime.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A Derivation of analytical solution of linear fluid flow in porous media
One dimensional fluid flow is considered in a sample with dimensions of L Â H Â W. Initial pore pressure of the sample is set equal to P i . Boundaries have constant pressure. Pressure at one end of the sample is P 1 and on the other end is P 2 . Equation 38 shows the initial condition. Equations 39 and 40 show boundary conditions: Equation 41 shows the linear form of pressure diffusion equation in porous media (Donnez 2012): a 2 is called hydraulic diffusivity and is equal to fluid mobility divided by fluid storability (Donnez 2012): Equation 41 can be solved by separation of variables to get an expression for pressure function that its variables are time and location. Equation 41 can be solved easily for homogenous boundary condition, i.e., boundary condition values set to zero. But in this situation one or both of boundary values can be a pressure above zero. To overcome this problem and convert it to a homogenous boundary condition equation, the pressure function can be assumed to be composed of two parts of one being time-independent and the other part being time dependent (Gonzalez-Velasco 1996): In Eq. 43, vðxÞ is the steady state pressure equation that is independent of time and aids to change the boundary conditions of the problem to homogenous boundary conditions. wðx; tÞ is the transient part of pressure function and its value will change with respect to time.
Because vðxÞ is part of the pressure function equation, it should satisfy the diffusivity equation as well. vðxÞ is independent of time and its derivative with respect to time is zero as shown in Eq. 44: Replacing Pðx; tÞ by vðxÞ in Eq. 41 and using Eq. 44 gives: The solution of Eq. 45 is: Rewriting the boundary conditions for vðxÞ; P 0; t ð Þ ¼ v 0 ð Þ ¼ P 1 and P L; t ð Þ ¼ v L ð Þ ¼ P 2 . Values of c 1 and c 2 can be determined by these two boundary values (Gonzalez-Velasco 1996) Steady state part of the pressure equation is determined. The next step is to find the transient part of the pressure equation. The boundary values for the transient part are calculated as follows (Gonzalez-Velasco 1996): As w x; t ð Þ is also part of the pressure equation, so it should satisfy the partial differential pressure diffusivity equation: This equation can be solved by the method of separation of variables. If w x; t ð Þ is denoted as a multiplication of two functions that one of them is only dependent on x and the other is dependent on time then (Schröder 2009): Partial derivatives of w x; t ð Þ with respect to time (t) and position (x) are: Substitution of Eqs. 52 and 53 into Eq. 50 gives: If either of XðxÞ or TðtÞ is zero, then the solution of wðx; tÞ is the trivial solution of w x; t ð Þ ¼ 0. So both of these functions are different from zero and both sides of Eq. 54 can be divided by X x ð ÞTðtÞ to give: Sides of Eq. 55 are independent of each other. The left hand side is a function of time (t) and the right hand side is a function of position (x). So in order for equation to hold for any t and x, both sides should be equal to a constant (Schröder 2009). If sides are equated to -k (-k can be any positive or negative number but the negative sign makes calculations simpler) it will result in: And, So Wðx; tÞ decomposed to a homogenous second order ordinary differential equation for XðxÞ and a homogenous first order ordinary differential equation for TðtÞ. Solving for boundary conditions: T t ð Þ ¼ 0 will make both of the conditions to satisfy but it also causes w x; t ð Þ equation to be zero which is a trivial solution and not the desired solution. Therefore; To solve Eq. 56 there are three options for k: Option 1: k < 0 Option 2: k = 0 Option 3: k > 0 Three options will be considered individually to decide which option will give the correct answer: Option 1: k < 0 k ¼ Àd 2 so the characteristic equation is (Kreyszig 2010): r 2 À d 2 ¼ 0 and roots are r ¼ AEd.
This gives the general solution as (Kreyszig 2010): Putting the first boundary condition from Eq. 58 into Eq. 60 gives: Putting the second boundary condition from Eq. 59 into Eq. 60 gives: ðe dL þ e ÀdL Þ 6 ¼ 0 and for Eq. 62 to hold, D 1 should be zero which causes D 2 to be zero as well. This will result in X x ð Þ ¼ 0 which in turn causes w x; t ð Þ ¼ 0 which is a trivial solution. So this means that option 1 is not valid and k cannot be a negative value.
Option 2: k = 0 This option will result in Eq. 63 for X x ð Þ: Applying the first boundary condition will results in: And applying the second boundary condition will result in: This will result in X x ð Þ ¼ 0 which in turn causes w x; t ð Þ ¼ 0 which is a trivial solution. So this means that option 2 is not valid and k cannot be zero and only option 3 is left.
Option 3: k > 0 k ¼ d 2 so the characteristic equation is (Kreyszig 2010): r 2 ? d 2 = 0 and roots are r ¼ AEdi. So based on the roots of the characteristic equation the general solution is (Kreyszig 2010): By imposing the boundary conditions: . . .; np L . So based on the above three options being analysed, the value of k should be greater than zero. The values of k that satisfy Eq. 64 are eigenvalues of this equation and are: For each of eigenvalues there will be an eigen function. The nth eigen function is: Inserting the value of k from Eq. 65 into Eq. 57 will give: Solving Eq. 67 will give T(t) for nth eigenvalue: T n t ð Þ ¼ C n e À a 2 n 2 p 2 L 2 t ; n ¼ 1; 2; 3; . . .

ð69Þ
The summation of all solutions for different eigenvalues is also a solution of w x; t ð Þ. As a result w x; t ð Þ can be shown as: C n e À a 2 n 2 p 2 L 2 t sin npx L ; n ¼ 1; 2; 3; . . .
Appendix B Derivation of analytical solution of radial fluid flow in porous media in a finite circular reservoir Radial fluid flow is considered for a sample with external radius (R e ) and wellbore radius ðR w Þ. Initial pore pressure is P i . Pressure at outer boundary ðR e Þ is kept constant at P i while wellbore pressure changes to ensure constant wellbore flow rate. Equations 75 and 76 show inner and outer boundary conditions (Ahmed and McKinney 2011). Equation 77 is the initial condition.
The radial form of pressure diffusion equation in porous media is shown in Eq. 78 (Dake 1978): To simplify (78) and boundary conditions, the following dimensionless parameters are defined: R D is dimensionless radius, t D is dimensionless time and DP D is delta dimensionless pressure. Inserting Eqs. 79, 80 and 81 into Eq. 78, will simplify it to Eq. 82: Boundary conditions and Initial condition will simplify to: Equation 82 can be solved in the same way as linear form of the equation has been solved by assuming that the equation is composed of two parts of steady state and unsteady state.
S R D ð Þ is the steady state part of the equation and U R D ; t D ð Þis the unsteady state part of the equation. Steady state part of the equation should satisfy initial and boundary conditions as well as pressure diffusion equation. As the steady state solution is independent of time, the right hand side of the diffusivity equation is equal to zero.
Solving Eq. (87) gives: Inserting boundary conditions from Eqs. (88) and (89) into Eq. (87) results in C 2 ¼ lnðR De Þ and C 1 ¼ À1 and therefore; The steady state part of the dimensionless pressure equation is solved. The unsteady state part needs to be solved. But first its boundary and initial conditions need to be determined. Using the second part of Eq. 83: Inserting Eqs. 84 and 89 into Eq. 86, gives: Inserting Eqs. 85 and 91 into Eq. 86, gives: Equations 92 and 93 are the boundary conditions for the unsteady state equation and Eq. 94 is its initial condition.
Same as what has been done to solve the transient part of the pressure equation for linear flow, in here the concept of separation of variables is used to solve U R D ; t D ð Þ. U R D ; t D ð Þ can be shown to be a multiplication of two functions with one being dependent on time and the other dependent on space or radius.
By applying boundary conditions from Eqs. 92 and 93 to Eq. 95: ; which is not of interest. Therefore; As U D R D ; t D ð Þ is part of dimensionless pressure equation, it should satisfy the dimensionless partial differential equation: Partial derivatives of U D R D ; t D ð Þ with respect to time and space are: Substitution of Eqs. 99 and 100 into Eq. 98 gives: If either of T t D ð Þ or X R D ð Þ is equal to zero then U D R D ; t D ð Þwill be equal to zero which is not of interest. Therefore, both T t D ð Þ and X R D ð Þ are different from zero. If both sides of Eq. 101 are divided by T t D ð Þ Â X R D ð Þ, it gives: Sides of Eq. 102 are independent of each other. The left hand side is dependent on location or radius and the right hand side is dependent on time. In order for Eq. 102 to hold for any time and location, both sides should be equal to a constant. Sides can be equated to Ç . Same arguments can be made about the sign of Ç as before for the linear pressure diffusion equation. In here, it can be shown that Ç should be a negative value Àk 2 for equation to hold. So equating both sides to Àk 2 will give: Solving Eq. 103 gives: In Eq. 104, A is a constant.
Equation 105 can be solved in maple as shown below: dslove diff ðdiffðX R ð Þ; RÞ; RÞ þ 1 R Á diffðX R ð Þ; RÞ þ k 2 Ã X R ð Þ ¼ 0 X R ð Þ ¼ C 1 Besse1J 0; kR ð Þþ C 2 Besse1Y 0; kR ð Þ So the answer to Eq. 105 is: J 0 in Eq. 106 is first kind of Bessel function of order zero and Y 0 is second kind of Bessel function of order zero (Polyanin and Manzhirov 2008). To determine constants of the equation, boundary conditions form Eqs. 92 and 93 are applied to Eq. 106: Solving Eq. 107 for C 2 gives: Replacing above equation in Eq. 108 gives: If C 1 ¼ 0, it implies that C 2 ¼ 0 and as a result X R ð Þ ¼ 0 which is not of interest. Therefore; Knowing R De , Eq. 109 can be solved to get the value of k. There are infinite numbers of k that will satisfy this equation and they are called eigenvalues of this equation. Eigenvalues are represented as k n . For every k n; there is a C 1n , C 2n and as a result there is a X n R ð Þ: As every X n R D ð Þ is a solution for Eq. 95, their summation is also a solution for Eq. 95. Inserting Eq. 104 and 111 into Eq. 95 gives: Applying the initial condition from Eqs. 94 to 112 gives: C n can be shown to be (Muskat 1982): Inserting Eq. 114 into Eq. 112 gives: pJ 2 0 k n R De ð Þ k n ðJ 2 1 k n ð Þ À J 2 0 k n R De ð ÞÞ ðY 1 k n ð ÞJ 0 k n R D ð Þ À J 1 k n ð ÞY 0 k n R D ð ÞÞe Àk 2 Equation 115 is the final solution of the un-steady state part of pressure equation. Inserting Eqs. 91 and 115 into Eq. 86 gives the final solution to pressure equation (Muskat 1982): Dimensionless pressure is defined as: Wellbore flow rate at any time is constant and is equal to steady state flow rate. Steady state flow rate is: Inserting Eq. 118 into Eq. 81 gives: Inserting Eq. 119 into Eq. 116 gives: pJ 2 0 k n R De ð Þ k n ðJ 2 1 k n ð Þ À J 2 0 k n R De ð ÞÞ ðY 1 k n ð ÞJ 0 k n R D ð Þ À J 1 k n ð ÞY 0 k n R D Þ ð Þe Àk 2 n t D Equation 120 is the dimensionless pressure distribution across a circular finite reservoir with constant outer boundary pressure and constant wellbore flow rate.