A Review of the Dynamic Modeling Approaches for Characterizing Fluid Flow in Naturally Fractured Reservoirs

Fluid flow in fractured media has been studied for decades and received considerable attention in the oil and gas industry because of the high productivity of naturally fractured reservoirs. Due to formation complexity and reservoir heterogeneity, characterizing fluid flow with an appropriate reservoir model presents a challenging task that differs relatively from homogeneous conventional reservoirs in many aspects of view, including geological, petrophysical, production, and economics. In most fractured reservoirs, fracture networks create complex pathways that affect hydrocarbon flow, well performance, hence reservoir characterization. A better and comprehensive understanding of the available reservoir modeling approaches is much needed to accurately characterize fluid flow behavior in NFRs. Therefore, in this paper, a perspective review of the available modeling approaches was presented for fluid flow characterization in naturally fractured medium. Modeling methods were evaluated in terms of their description, application, advantages, and disadvantages. This study has also included the applications of these reservoir models in fluid flow characterizing studies and governing equations for fluid flow. Dual continuum models were proved to be better than single continuum models in the presence of large scale fractures. In comparison, discrete models were more appropriate for reservoirs that contain a smaller number of fractures. However, hybrid modeling was the best method to provide accurate and scalable fluid flow modeling. It is our understanding that this paper will bridge the gap between the fundamental understanding and application of NFRs modeling approaches and serve as a useful reference for engineers and researchers for present and future applications.


Introduction
Naturally fractured reservoirs (NFRs) can be defined as oil and gas reservoirs whose rock matrices have pores, pore-throats, and natural fractures, where fractures play a significant role in fluid flow and hydrocarbon recovery. NFRs can be found worldwide, containing substantial amounts of the world's hydrocarbon proven reserves [1][2][3][4][5][6][7][8][9]. In a study that has been conducted in 2019 by the International Energy Agency (IEA), British Petroleum and Schlumberger illustrated that around 60% of the world's oil proven reserves and 40% of the world's gas reserves are being trapped in naturally fractured reservoirs. Two different systems typically model these reservoirs; primary matrix and secondary fracture network systems [10][11][12][13][14][15][16]. The primary storage of hydrocarbons is provided by the rock matrix, which is characterized by low permeability and high storage volume. In contrast, the conductive path for fluid flow in porous media is governed by fractures that have high permeability but hold a small amount of fluid [17][18][19]. Fluid flow in fractures differs from that in a matrix. Most of the fluid in the matrix blocks moves relatively slowly, whereas, in fractures, fluid flows comparatively fast because fractures form high permeability paths [20].
Due to the decrease in hydrocarbon conventional resources with time and the high productivity potential of unconventional resources such as NFRs, many researchers began to find efficient means to model fluid flow in fractured porous media [21][22][23][24][25][26][27][28][29]. However, modeling and simulating multiphase flow and obtaining high recovery factors often depend on selecting the appropriate reservoir modeling process, for example, accurate and efficient modeling of a complicated fracture system and simulating its matrixfracture fluid exchange [30,31]. In NFRs, the hydrocarbon interflow within two different mediums, i.e., matrix and fracture, makes the modeling process very challenging [32]. These reservoirs are considered as one of the most common and crucial geological structures in oil fields [28,33,34]. Due to the high heterogeneity of their porosities and permeabilities, modeling and characterization of fluid flow are becoming more difficult throughout the entire world for both green and brownfields. Moreover, recent oil and gas field discoveries with Natural Fracture Networks (NFNs) are considered to be unique in terms of their geology and reservoir engineering issues, such as complications in Exploration and Production (E&P) operations [21,35,36].
In the reservoir engineering field, numerical simulation techniques are widely used to study flow and transport operations in subsurface fractured formations. The research about this topic has been advanced significantly in the past several decades due to its production potential [5]. There are various fluid flow characterization techniques; the choice of the method should be made according to the type and degree of heterogeneity of the reservoir as well as the primary objectives of the reservoir simulation study. However, the main remaining challenge in characterizing fluid flow in NFRs is selecting the appropriate reservoir modeling approach. The selection process is closely related to specific modeling problems and limitations discussed in this paper. Applying the suitable simulation methods can lead to accurate fluid production forecasts, optimized well-completion designs, the identification of the needs for an artificial lift for future processes, and making decisions about improved or enhanced oil recovery processes. This paper will cover several studies from a different perspective of previous work, including; reviewing the past and current naturally fractured reservoir modeling techniques and methods. It will also cover the applications of these models in characterizing fluid flow through such reservoirs. Furthermore, a comparison of fluid flow through fracture and matrix systems using different reservoir models will be discussed. In addition, governing equations for fluid flow will be presented. This paper will contribute to a better understanding of the simulation models available to characterize the process of fluid flow in fractured reservoirs. It will also highlight the role played by NFRs in reservoir production, recovery management, and optimization.

Types and Principles of Dynamic Modeling Approaches
Simulation modeling approaches can characterize reservoir depletion processes for a sophisticated, heterogeneous, and naturally fractured reservoir [37]. Most of NFRs are composed of different systems of matrices, fractures, and vugs. Thus, they have contrast in some properties, including porosity, permeability, and fluid transport functions. Therefore, identifying an accurate and efficient modeling technique to provide a detailed fluid flow characterization in fractured reservoirs presents a fundamental challenge [38][39][40][41][42]. Besides, fluid flow through fracture networks is a complicated inverse structure problem as multi-patterns, high-complexity of fracture distribution, and the solution's multiplicity defects are present [43]. As a result, a vital task is to establish the correct model to evaluate these reservoirs [44].

Single Continuum Fracture Model (SCFM)
To some extent, a single continuum fracture model, i.e., Single Porosity Model (SPM) or Single Permeability Model (SKM), can represent fracture systems' complexity. SCFM can be applied when the fracture pattern is scale-independent, i.e., when fractures have no orientation and the matrix blocks are small. It can also be used when the productivity is mainly dominated by flow movement rather than matrix fracture interaction. The single continuum modeling approach is also used to characterize fractured reservoirs with united fracture and matrix properties as well as reservoirs with unique and effective permeability. In this type of modeling, matrix-fracture flow is often considered a weak point treated through transfer functions. The most proposed solution to simulate NFRs using this method is based on Warren and Root Pseudo-Steady-State (PSS) model and its modifications [51,52]. In a reservoir with a defined boundary (closed reservoir), PSS refers to the flow region after the pressure wave had propagated through all the reservoir boundaries.
In this approach, fracture sets might be represented by the change in fractured medium permeability. Thus, the magnitude and the tensor direction of the system permeability may vary depending on fracture network properties [53]. Furthermore, in single continuum modeling, the reservoir properties might change from point-to-point and might be divided into a number of representative volumes and bulk macroscopic values. Using a single continuum fracture model requires averaging some of the fracture properties over the volume, often known as block-based permeability tensors [54][55][56][57]. As a consequence of this averaging process, some properties of separate fractures might be lost. In the numerical simulation of fluid flow in NFRs, each grid block contains two primary fracture and matrix variables. Therefore, capturing all fracture systems' complexity using single continuum modeling is highly difficult [18,35].
The single continuum approach is desirable in some particular cases because of its reduced level of computational requirements and model complexity. Furthermore, in these models, the simulation run times are very short, with no additional terms in comparison with a non-fractured reservoir model. Hence, the fluid flow exchange functions and simulation techniques are straightforward. However, some disadvantages of using this technique may include difficulties in providing a detailed simulation model for the complex transport processes in NFRs due to restrictions caused by the simplifications of SCFM. Furthermore, in some simple strategies, such as the pressure evolution process in tracer transport or compressible flow, the fluid flow dynamic behavior is not adequately represented by averaging the flow alone [18,35,53].

Dual Continuum Fracture Model (DCFM)
Dual continuum models may include dual-porosity (DPM) or dual permeability (DKM) models. Barenblatt et al. [58] introduced the continuum modeling concept. Then it was first applied by Warren et al. [59], where fractured reservoirs were presented with two independent and overlapping continua. The first medium is called the matrix, which denotes to fluid flow within the deformed element, i.e., the matrix domain, whereas fluid flow was simulated by the other component, i.e., fracture. Later, different studies extended these models to simulate multiphase flow [60,61]. Among a number of different studies using various modeling approaches for fluid flow characterization [62][63][64][65], it was stated that dual continuum modeling approaches are well suited to quantify hydrocarbon flow in NFRs.
A dual-continuum system consists of two continua; the globally permeable mobile continuum and the immobile continuum with negligible global permeability, which transports fluid between itself to the mobile continuum [66][67][68][69]. A dual continuum model's behavior might be characterized by the mass transfer rate of the mobile and immobile continua [67]. In addition, DCFMs are used only for matrix and fracture systems where oil is stored in matrix blocks and transmitted through fractured paths, for instance, in type II reservoirs. The overall conductivity of type II reservoirs is provided by fractures, whereas the rock matrix supplies the overall storativity. In type II reservoirs, the matrix has sufficient permeability to provide conductivity for fluid flow from the matrix to fractures. The initial rate for this type of reservoir can be high even in low matrix permeability conditions. When the communication between matrix and fractures is poor, type II reservoirs might be subjected to difficulties during secondary recovery mechanisms [70].
In dual continuum modeling, matrix and fractures are dealt with as two separate continuous mediums. These models can also be used when the fracture connectivity is high and when the matrix-fracture interflow is significant [71]. The advantage of the dual continuum models vs. other approaches is that DCFMs can be extended to multiphase flow without considering the fracture orientation complexity. The dual continuum models, such as dual-porosity models, assume a uniform distribution of the fractures throughout the reservoir [72]. Moreover, with dual models, the matrix-fracture interaction can be achieved in NFRs. However, this may lead to a longer running time than single continuum modeling approaches.

Dual Porosity Model (DPM)
According to Moinfar et al. [13], dual-porosity models involve the discretization of fractured reservoirs into two domains, i.e., matrix and fracture. Therefore, all reservoir grids contain fracture, matrix pressure and saturation properties. In this modeling technique, the interconnected fractures provide the fluid flow paths to producers and injectors. In contrast, the low permeability and high storability matrix act as a fluid storage sink. In DPM, within a gridblock, each fracture cell is connected to its corresponding matrix cell through a single exchange term. However, there is no direct communication between neighboring matrix blocks, where they are linked only through fractures. Dual-porosity models can simulate gravity drainage effects by modifying the source term to account for gravity. However, if the fracture and matrix depths are equivalent, these models are unable to model the impact of the gravity drainage [73,74] Double porosity models have been widely implemented in many full-field scale reservoir simulation studies [48]. To model fluid flow in large-scale fractures, several studies by Bourbiaux et al. [75], Gilman [76], Pruess [77], Duguid et al. [78] presented various matrix fracture transfer functions. The dual porosity approach was also used to provide an analytical and numerical solution in some studies [79,80]. In the dual-porosity modeling approach, a Transfer Function (TF) is used to exchange the hydrocarbon flow between matrix and fracture systems. In contrast, a shape factor is used to describe the fracturematrix surface area [81]. Dual-porosity models might approximate the effect of the gravity mechanism. DPM uses the assumption of segregated fluid, which is similar to the vertical equilibrium option in many single-porosity simulators.

Dual Permeability Model (DKM)
Dual permeability models are introduced when there is a need to simulate the matrix continuity in fractured reservoirs. Furthermore, in dual permeability models, fluid flow is allowed to occur between matrix grid blocks contributing to more realistic and accurate fluid flow modeling than in a dual-porosity model. In addition, dual permeability models have been utilized in many commercial reservoir modeling simulators [82]. DKM allows the use of production mechanisms such as gravity and capillarity mechanisms as well as capillary continuity or discontinuity of matrix blocks [83].

Dual Porosity/Dual Permeability Model (DPDP)
DPDP models can be used in combination to characterize fluid flow processes in NFRs with more comprehensive and accurate characterization results. DPDP models are the most widely used idealization to simulate large-scale petroleum reservoirs. In addition, these models have been used in the industry to history match numerous fractured reservoirs and predict their future performance [84]. To effectively model these processes, the capillary, viscous, gravitational effects must be taken into consideration. The fractured reservoir's type (geometry) and its recovery mechanism define the significance and interaction between these forces. Viscous displacement from the matrix in these models is generally neglected but, Gilman et al. [85] presented a method to account for this displacement mechanism.
In DPDP simulators, fractures and matrix gridblocks are regarded as two solid mediums linked to one another. Thus, flow and medium parameters are determined in each of the mathematical grids. Three different paths can characterize the fluid flow in these models; fluid flows through fractures, the matrix (the dualpermeability model), and finally, between matrix blocks and fractures. The fluid flow along the fractures is defined by filtration equations, particularly by black oil simulators. The additional sink or source type member is used to characterize fluid flow between matrix and fracture systems.
In dual-porosity dual-permeability models, the fracture cells typically have much smaller porosities and smaller pore volumes than the matrix cells, even though both continua cells have the same bulk volume. Smaller fracture cell size coupled with the higher fracture permeability usually produces different flow time-scales between the two continua. It makes the numerical simulation process of NFRs challenging and more computationally expensive, mainly due to the smaller time-step requirements. It is crucial to separate the fractured medium's absolute and effective permeabilities [45]. For full-field scale flow simulation studies, the use of dual-porosity/dual-permeability approaches requires averaging flow behavior over a large volume, resulting in some modeling limitations to represent the fluid flow in fractures and from the matrix to fractures [84].
Previous research about fluid flow using dual-porosity/dual permeability models was only limited to simplified orthogonal natural fractures or 2D models [86] and concentrating mainly on the hydraulic problem. In addition, the effect of natural fractures with different dip and strike angles on the hydromechanical behavior of a fractured reservoir is insufficient to be characterized under the assumptions of 2D models. Therefore, DPDP models have been further enhanced (EDPDP) by Rueda et al. [87] to represent a fractured porous formation more realistically using a 3D hydro-mechanical formulation. This modeling approach incorporated multi-block domains through multiple fracture sets with various scales, arbitrary orientations, permeabilities, spacing, and sizes. It also included fractures that can be penetrated by secondary and tertiary fractures. The EDPDP model was implemented in an in-house framework GeMA (Geo Modelling Analysis) using the finite element method [88]. These models were then compared to the discrete fracture method to assess their accuracy. The results suggested that the proposed model was accurate enough to model full-field scale highly fractured rocks.
Since dual-porosity continuum models assume uniform matrix and fracture properties throughout the reservoir, they may not capture the fractured reservoirs' actual complexity [89,90]. To overcome the limitations of the dual-continuum approaches, Ci-qun [91], Ci-qun [92] proposed an improvement to the issue by considering two matrix systems with different properties or by considering two fracture systems with different properties in addition to the matrix. The latter component is normally referred to as the dual fracture model. Ci-qun [91], Ci-qun [92] model was used for the radial flow of slightly compressible fluids through a triple-porosity reservoir under pseudosteady state inter-porosity flow. Another tripleporosity model was introduced by [40], where two different geometrical configurations were used (strata model and uniformly distributed blocks model). This model has included two matrix systems and a single fracture under unsteady state inter-porosity flow. Reference [93] investigated the transition zone behavior of the radial triple porosity system. They extended the [40] strata (layered) model by allowing the matrix systems to have different properties and thickness. In addition, pseudosteady and interporosity unsteady flow conditions were taken into account.
Triple governing equations are used to describe fluid flow and transport processes in triple triplecontinuum approaches separately. The obtained conceptualization results are presented in the same form as that for a single porous medium, where a set of partial differential equations are solved for each continuum [94]. Pruess [95] stated that the triple-continuum models are not restricted to the orthogonal idealization of the fracture systems. These triple-continuum models can handle the complexity and heterogeneity in the matrix and fracture systems [96]. A radial triple-continuum model was developed by Liu et al. [39] that contained fractures, matrix, and cavity medium. Their solution was an extension of Warren et al. [59] solution. Unlike previous triple-porosity models, fluid is allowed to flow from the matrix and cavity systems under pseudosteady state to the fractures and then to the producing well; thus it is called triple-continuum. To more accurately characterize the geological structure of the coal seam and its production forecasts, a triple-porosity/dual-permeability model (TPDP model) was developed. In the TPDP model, the matrix consists of interconnected macro-pores and discontinuous micro-pores. Thus, there are three pore systems: the micro-pore system, the macro-pore system, and the fracture system. Apart from the micro-pore system, the other two are continuous [97].
To capture the flow mechanism in NFRs, several multi-porosity models have been developed and used. Al-Ghamdi et al. [98] introduced two sub-triple-porosity models with a single matrix and two fracture systems, where macro-fracture was more permeable than micro-fracture. The macro-fractures and/or micro-fractures are allowed to flow to the well under pseudosteady state flow. Dreier [99] enhanced the triple-porosity dual fracture model originally presented by Al-Ghamdi et al. [98]. The flow between micro and macro fractures was represented by a transient flow condition, whereas the matrix and micro-fractures flow were still under pseudosteady state conditions. Reference [94] characterized the flow and transport of tracers and nuclear waste in the unsaturated zone of Yucca Mountain using a triple-continuum modeling approach. Large fractures, small fractures and matrix blocks were all present in the system. This model was useful in estimating reservoir parameters. The method was validated by comparing the analytical solution with a numerical simulation model.
A more recent study introduced a triple-porosity model [100] to effectively model and capture rock matrix and vugs and the microfractures connecting them. This research was conducted to understand the fluid flow behavior of fractured carbonate reservoirs and capture their heterogeneity. In this study, a Discrete Fracture and Vug Model (DFVM) was applied to model macrofractures and macrovugs. Whereas microfractures and microvugs were modelled with the triple-continuum concept. Macrofractures and the boundaries of macrovugs were discretized into several elements. The Boundary-Element Method (BEM) was used to handle flow into macrofractures and macrovugs. The flow within macrovugs was assumed to be pseudo steady-state to decrease the computational burden. The model was then validated with an analytical solution.

Discrete Fracture Network Model (DFNM)
The fracture network that disregards the host medium (matrix system) is commonly defined as a discrete fracture network. This modeling approach can be used when the entire storage and conductive paths of fluid flow in the reservoir are governed by fractures that can be represented explicitly. Therefore, fluid flow is only restricted to fractures neglecting the flow-through matrix blocks and isolated fractures. DFNMs can also be used as a generalization model to characterize fractures in heterogeneous and low-permeable reservoirs [47,53,72]. This approach establishes an interconnected fracture network surface by mapping the fracture planes spatially on three-dimensional grids. Gilman [84] reported that a reservoir rock volume could be regarded as a matrix block when it is three-dimensional and bounded by fracture planes.
The discrete fracture network model considers all fracture properties to model an explicit, realistic, and accurate fracture geometry and complex flow in a fracture network [101][102][103][104][105][106]. Explicit approaches have found merit in a number of engineering applications, including; fault reactivation [107,108], hydraulic fracture propagation [109], the interaction between hydraulic and natural fractures [110][111][112], fluid flow in naturally fractured reservoirs [113], among others. DFNM approaches all rely on matching fracture and matrix grids in the sense that a fracture element coincides geometrically with co-dimension-one mesh entities, i.e., faces of matrix grid elements [106].
Constructing a DFNM has been widely presented in previous literature, e.g., [114][115][116][117][118][119]. For each fracture system, it is crucial to input several fracture parameters such as the independent fracture location, the fracture cluster, length, aperture, orientation (angle of dip and dip azimuth), permeability, and volumetric fracture intensity. These fracture parameters are imported directly into the model without any averaging. The accuracy of these models depends on the reliability of the data and its source. Fractures in DFNMs are dealt with as 3D objects. These models are highly affected by fracture orientation length, intensity, and transmissivity [120]. Fracture connectivity is crucial for assessing the behavior of hydrocarbon flow in NFRs [121].
A discrete fracture network might be derived from different sources and scales ranging from kilometers (e.g., surface seismic data) to centimeters (e.g., borehole image and core analysis) [122]. Some methods may derive discrete fracture networks using a statistical approach, such as defining fracture orientation and intensity using wide azimuthal seismic data [123,124] or using 3D vertical seismic profiling [125]. Other methods generate explicit discrete networks where the fracture networks' statistical information and the precise details of each fracture patch, such as location, shape, and size are presented [126].
The advantages of using discrete fracture network models are; they allow the quantification of some of the flow and transport phenomena, which may not be captured using continuum modeling techniques. Discrete modeling can also take into account the impact of individual fractures on fluid flow and solute transport. The fracture geometry of NFRs can be realistically represented when a DFNM is directly used to model the production layer connectivity through fractures. Furthermore, DFNM has a great potential to resolve well-completion geometry issues to improve hydrocarbon production and recovery. Also, breakthrough times and spatially distributed fractures can be captured more accurately. DFNMs can also be used to estimate fracture behavior and connectivity based on well test analysis and single or multiple well interference tests. Moreover, any single fracture is represented discretely rather than presenting them as a group of regularly spaced fractures within the matrix cubes [18,35,47,84].
On the other hand, some drawbacks of DFNMs have made them inapplicable in complicated reservoir modeling. If the modeling process involves many fracture systems, a large number of grids are required, which results in extensive computational resources and running times to solve the system equations compared with other approaches. This method is also inappropriate to be used in multiphase flow modeling and full-field scale simulation studies. Since discrete fracture network modeling only considers fluid flow in fractures, a modification should be applied to forecast the permeable matrix's fluid flow behavior. Finally, reservoir zones far from the wellbore lack reliable information, making DFNM very uncertain. Moreover, the inter-porosity flow adds more complexity, requiring extensive computational power and scale processes that have not been developed for modeling multiphase fluid flow [18,84].

Unstructured Discrete-Fracture Model (USDFM)
This approach was initially presented by Karimi-Fard et al. [105]. This method's gridding system is unstructured, where a lower-dimensional technique is applied to grid the DFNM models. Rock matrix modeling is based on 3D polyhedral cells, whereas fractures are characterized by a subset of the 2D interfaces separating the grid blocks. Some knowledge about the neighboring control volumes (a connectivity list) and the transmissibility at each connecting point is required to account for fluid interflow between adjacent grids [127,128]. Reference [113] presented a higher-order numerical model for fractured media (2D and 3D unstructured gridding) based on an efficient algorithm for two-phase compositional flow. This approach accommodated all types of finite elements, including quadrangular and triangular elements in 2D, and hexahedra, prisms, and tetrahedra elements in 3D. Their algorithm allowed for flow simulation in fractures considering all range of permeability values as opposed to other models where low permeable fractures affected the accuracy of the results.

Embedded Discrete Fracture Model (EDFM)
Gridding fracture objects' challenges and limitations are entirely circumvented using the EDFM technique based on discrete fracture conceptual modeling, where fluid flow in both matrix and fracture systems is discretized separately [129]. Fluid flow through the matrix and fractures within the matrix cell is proportional to the pressure difference between the two mediums, where fracture geometry is reflected by adding a constant. In this approach, the discrete model is designed using a conventional structured grid. Additional control volume of fractures is created by simulating the fracture intersection with the matrix grid cells. Furthermore, matrix-fracture flow can be achieved by transport parameters between fracture networks and the matrix blocks [130]. A significant advantage of EDFMs is their "non-intrusive" nature, allowing them to be easily incorporated into reservoir simulators that support non-neighbor connections (NNCs) [13,131]. The NNCs occur once all the fractures are connected with the rest of the reservoir. Several non-neighbor connections can be used in the embedded discrete fracture modeling, including; the connection of two intersecting fracture segments from different or same fractures and fracture segment connection with its parent matrix grid blocks [132].
An efficient approach was presented by [133] to simulate natural fractures of complex geometries and non-planar hydraulic fractures. After which, it was compared with local-grid-refinement (LGR) models and a semi-analytical solution. NNCs were used with EDFM. In this research, discrete fractures were installed into the computational domain using the non-intrusive feature of EDFM and the simulator's original functionalities. Reference [134] extended the approach to account for complex fracture network simulation using corner-point grids (CPGs). A general geometrical algorithm was introduced to find the matrix-fractures intersection points. Moreover, the transmissibility factor was further enhanced for the connections and intersections between fractures. Reference [135], developed a new approach to study the dynamic behavior of complex fractures by reformatting the EDFM used in previous models and linking the gridblock permeabilities to EDFM transmissibility factors. Therefore, by varying the permeability of fracture gridblocks with time, complex fractures' dynamic behavior can be conveniently simulated. This method has no difficulty in gridding around fractures because the fractures and matrix are separately gridded. In addition, it is more convenient to simulate large number of fracture networks or complex fracture geometries.
The development of naturally fractured reservoirs, especially shale gas and tight oil reservoirs, exploded in recent years due to advanced drilling and fracturing techniques. However, complex fracture geometries such as irregular fracture networks and non-planar fractures are often generated, particularly in natural fractures. Accurate modeling of production performance from reservoirs with such geometries is challenging. Therefore, it is essential to use EDFM to effectively define the complex fracture geometries and identify their properties [136]. Embedded discrete fracture models are proved to possess high computational efficiency and flexibility to model complex fractures. Reference [137] introduced a selection of field-scale applications of these models to characterize shale reservoirs accurately. The results indicated that the EDFMs are practical to be used in conjunction with history matching and sensitivity studies that include complex fractures.

Geological Discrete Fracture Network Model
Fractures in this approach are mapped from the rock outcrops exposure. These models have several roles in the oil and gas industry, including understanding the process of fracture formation, tectonic stresses interpretation, and derive the statistics and scaling of fracture distribution [138,139]. These models can also provide more details about fracture connectivity and multiphase fluid flow [140]. Detailed field mapping may determine fracture apertures, which are often considered constant or follow statistical distribution [141]. The geological DFNM technique can preserve the feature of natural fracture (e.g., curvature and segmentation). Besides, the characterization of complex fracture networks can be done in an unbiased way. Nonetheless, this method is restricted to 2D models and may not be applicable for 3D analysis involving dipping fractures.

Stochastic Discrete Fracture Network Model
Stochastic modeling approaches use statistics from limited sampling intervals to perform 3D comprehensive measurement of the natural fracture system. Stochastic simulation of fracture systems is an accurate approach to build a set of dense and complex fracture networks. In 2D, fractures are assumed to be straight lines, whereas, in 3D, they are assumed planar discs or polygons. Other properties such as position, orientation, and aperture are treated as independent random variables following specific probability distributions [142]. Furthermore, in this approach, fractures' location is also assumed to be random, where corresponding probability density functions are used to sample the fracture geometrical attributes [143]. Fracture frequency might be related to fracture density or fracture intensity. While fracture density is a measure of fracture numbers, fracture intensity measures the total fracture persistence. The fracture spacing is described in terms of fracture frequency [144]. Fracture networks are grouped based on their orientation and characterized as uniform, normal, or fisher fractures [145].
Stochastic DFNMs are associated with significant uncertainties because of several assumptions, including; homogeneous spatial distribution, neglecting the correlations of different geometrical attributes, and simplifying fracture shape [140]. This can be improved by (1) spatial clustering process where inhomogeneity of fracture distribution is considered, (2) the correlation between fracture attributes, and (3) characterization of fracture intersections. It is more suitable for highly persistent fracture systems and accurate detection of fracture locations in the stochastic generation to use the spacing distribution approach [146,147].

Geomechanical Discrete Fracture Network Model
Geomechanical DFNMs are mainly dependent on geometrical stress and deformation to simulate fracture network evolution and represent fracture growth physics. Several numerical methods have been presented in this field, and the most frequently adopted method is the one based on Linear Elastic Fracture Mechanics (LEFM) [138,139]. Fracture patterns in the LEFM can be modeled iteratively by four main stages; (1) initial flaws development to imitate the process of natural fracture created from the microcracks, (2) perturbed stress field calculation in the rock, (3) stress intensity factor (KI) derivation at the fracture tip, and (4) fracture propagation which satisfies a growth criterion [148,149]. These models have the merit of connecting fractures geometry and topology with the conditions and physics of their formation. Moreover, Geomechanical DFNMs can automatically correlate each fracture's geometrical properties by the governing physics. However, this concept is associated with some degree of uncertainty and difficulty developing a consistent fracture pattern with the actual system [150]. A general numerical model for discrete fracture networks (DFN) is usually employed to handle the observed wide variety of fracture properties and the lack of direct fracture visualization. These models generally use fracture properties' stochastic distribution based on sparse and seismic data without any physical model constraint [151]. The usual approach in reservoir numerical modeling for fractured reservoirs is to use the statistical distribution of different fracture network parameters and calibrating them to local well log, seismic, and outcrop data [152]. Apart from the LEFM approach, fractured reservoirs can be simulated using other numerical techniques. Cowie et al. [153], Cowie et al. [154] developed a lattice-based rupture model to simulate the anti-plane shear deformation of a tectonic plate and the spatiotemporal evolution of a multifractal fault system. Tang et al. [155] used a damage mechanics FEM model to simulate the evolution of parallel, laddering, or polygonal fracture patterns formed under different boundary conditions, i.e., uniaxial, anisotropic and isotropic tectonic stretch, respectively. Spence et al. [156] employed the discrete element method (DEM) to simulate the fracture pattern development in a sedimentary sequence embedded with stratified nodular chert rhythmites. Asahina et al. [157] coupled the finite volume multiphase flow simulator (TOUGH2) and a lattice-based elasticity and fracture model (Rigid-Body-Spring Network) to simulate the desiccation cracking in a mining waste material under a hydromechanical process. The finite-discrete element method (FEMDEM) has also been used to simulate the formation and evolution of geological fracture patterns [158,159].

Hybrid Fracture Modeling (HFM)
Due to some limitations in applying each modeling approach individually, combining the different modeling methods to characterize fluid flow in fractured reservoirs is beneficial. For instance, integrating continuum modeling with DFNM is undoubtedly helpful in providing more accurate characterization results. Integrated models can accommodate various modeling scales and multiple data sources [18,31,111,160,161]. Tab. 1 summarizes the available modeling approaches in terms of description, area of application, advantages, and limitations.
The fractured/vuggy reservoir contains a low permeability rock matrix with numerous multiscale nonhomogeneous vugs and fractures of various diameter sizes ranging from meters to millimeters. Hybrid models are used to simulate multi-scale fractures and vugs by integrating continuum and discrete-fracture model concepts. Microfractures and vugs are represented with double-and multi-porosity concepts, whereas micro-fractures and vugs are simulated with DFNMs. In hybrid models, fine grids can describe the areal heterogeneity and the shape of macro-vugs. However, these models are associated with a heavy computational burden [162][163][164].

Description
• United and averaged properties.
• Reservoirs with single effective permeability.
• Matrix and fracture are regarded as two separate continuous mediums.
• All fracture properties are considered without any averaging.

Area of Implementation
• No fracture orientation.
• Matrix block is small. • The productivity is mainly dominated by flow movement rather than by matrix fracture interaction.
• The fluid is stored in the matrix and transmitted through the fracture • Fracture connectivity is high enough.
• Resolve well-completion geometry issues.
• Estimation of fracture behavior and connectivity.
• Reservoir porosity and permeability are due to fractures.
• Short run times.
• Establishment of the matrix-fracture interaction.
• Model extension to multiphase flow without fracture orientation complexity.
• Allow the quantification of flow and transport phenomena.
• Consider the impact of individual fractures.
• Breakthrough times and fracture geometry can be captured more accurately and realistically.

Limitations
• Difficulties in providing detailed modeling of the complex transport processes.
• Longer simulation running times than in single continuum models.
• Extensive computational resources and running times.
• Not applicable for multiphase, full-field scale simulation studies.
• Deep reservoir zones lack reliable information, which increases uncertainty.
• Complexity due to the inter-porosity flow.
The hybrid modeling technique can forecast short-term and long-term production scenarios (history matching) in NFRs. For instance, various fractured gas reservoirs with different forecasting horizons were modeled and matched by integrating multiple techniques, including; decline curve analysis, physics-based principles, and deep learning models. The outcomes suggested that combining more than one approach is highly beneficial for NFRs production prediction. However, HFM usually depends on production pattern complexity, data availability and accuracy, and the forecasting horizon. Thus, the accuracy of any reservoir model is a crucial factor in identifying the performance of the hybrid modeling technique [165]. A coupled simulation technique under isotropic and anisotropic in-situ stress conditions was introduced by [166] to study fracture propagation in a fractured formation. The approach combined the Finite Element Method (FEM) and Discrete-Element Method (DEM) using a reservoir simulation software known as TOUGH2.

Early Development of Fluid Flow Modeling Concepts
To analyze fluid flow in a slightly compressible, homogeneous, and double-porosity reservoir Barenblatt et al. [58] introduced a system of equations for fracture and matrix systems and a transfer function between them. A pseudo-steady state condition was used, neglecting the inertial effect. However, solving these equations was inconvenient and complicated in fluid flow characterization, raising the attention to develop a more realistic model to characterize fluid flow in NFRs. One of the classical applications of the dual-porosity technique was studied by Warren et al. [59]. An idealized model, as shown in Fig. 2, has been developed by Warren et al. to analyze the characteristic behavior of a permeable porous medium (a naturally fractured or vugular reservoir) that contained regions with a massive contribution to the pore volume of the system, but a minor contribution to the flow capacity. Using unsteady state (USS) flow and assuming a homogeneous reservoir, the well testing performance for this reservoir model has been investigated analytically where a sugar cube representation was used [59]. USS refers to flow parameters such as velocity, pressure, and density for each point that are time-dependent.
Odeh [167] Figure 2: The idealization of the heterogeneous naturally fractured reservoir, modified from [59] and Root model) may not be applied to every fractured reservoir. Kazemi [168] presented an ideal theoretical model to describe flow in a uniformly distributed fractured medium, as shown in Fig. 3. His study was conducted to modify and develop Warren and Root's model. Furthermore, in Kazemi's model, the reservoir was layered, finite, and circular, which contained two distinct porous regions (matrix and fracture) and a centrally located well. Fracture possessed low fluid storage but high flow capacity, whereas matrix blocks had high fluid storage and low flow capacity. He found an approximate value for matrix permeability by combining interference and build-up test results. In another study conducted by Reiss [169] to model fluid flow in a fractured porous medium, the reservoir was modeled as a set of vertical match stick columns, as shown in Fig. 4 [170].  Since the 1980s, discrete fracture network models were used to investigate fluid flow behavior in fractured reservoirs. The early fractured reservoir models were deterministic and straightforward, where fractures were fixed and predefined as a group. More complicated and stochastic models were soon proposed to characterize hydrocarbon flow through fractured reservoirs accurately. Many authors have presented different discrete fracture models. One such model is the Dershowitz polygon fracture conceptual model, Fig. 5. Polygonal fractures were defined on fracture planes by the process of Poisson lines resulting from the intersection of fracture planes. Fracture closing at fracture intersections can be accurately modeled using this approach. However, this model requires an indirect determination of fracture size and closing through the fracture planes' intensity [143,171]. Tab. 2 indicates a summary of the early main modeling concepts.

Dynamic Model Description
The dynamic simulation model can be defined as the reservoir's mathematical model developed to predict reservoir performance using specific computer software. These computer programs use sophisticated numerical methods such as the finite volume method, finite difference method, and finite element method to solve flow equations occurring during the hydrocarbon production phase. Usually, the dynamic model construction is an iterative process based on calibrating fracture properties against real data that provides information on fractured reservoirs' fluid flow behavior. The most important feature of a dynamic model is that it combines static model, pressure, saturation properties, well locations, and geometries with dynamic properties (production and well test data) to obtain production profiles vs. time. Furthermore, dynamic flow analysis can couple fracture properties with time, pressure, and flow to characterize fluid flow accurately in fractured systems. Hence it helps in selecting the appropriate development strategy to predict the best reservoir performance. The primary sources for dynamic data are well tests, production logging tests, and well performance information. A generated fracture network might also be an input for a flow simulation model [35].
A dynamic model can be used several times to simulate the entire life of a reservoir based on utilization schemes and operating conditions to optimize its depletion. Therefore, improved characterization techniques are required to accurately model heterogeneous fractured reservoirs and predict their fluid flow characteristics in regions distant from producing wells. Realistic models such as DFNM might lead to identifying high productivity zones (sweet spots), causing a positive effect on the reservoir production strategy and economic evaluation [3,172].
The dynamic permeability is achieved by the fracture's deformation (opening or closure), which is caused by the change of stress applied on the fracture planes due to oil and gas production or water injection [173]. The dynamic fracture growth is typically considered to be relatively slow. In a low permeability reservoir, it can run through the whole field development process. Dynamic fractures play a vital role in fluid flow modeling. In the optimistic case, they can highly improve the permeability of the reservoir and its production recovery. Still, in other cases, they might increase the reservoir heterogeneity leading to water channeling and strong water flooding [48].

A Comparison of Fluid Flow through Fracture and Matrix Systems Using Reservoir Modeling Approaches
Primary fluid flow paths in a typical fractured reservoir are mainly provided by fracture networks owning to their high permeability values. In contrast, the reservoir's main storativity is granted by the porous matrix volume due to its low permeability value [36]. In a fractured medium, natural fractures play a crucial role in controlling hydrocarbon flow. In any reservoir simulation model, it is critical to explore the static properties such as hydrocarbon in place and the dynamic behavior (e.g., oil production rate) realistically by understanding the fluid flow paths in different models [174]. Fundamentally, each reservoir modeling approach can be distinguished based on the capabilities of flow and storage for both fractured mediums, i.e., matrix and fracture [47]. In NFRs, each phase of fluid flow is distributed unevenly and governed by complex processes such as phase isolation [175]. A preferred fluid flow path is created when the fracture aperture and length become higher than a specific value. Fracture aperture and permeability are not related linearly but might be correlated between fracture transmissibility and flow. Therefore, when fractures dominate the fluid flow, the matrix-fracture permeability ratio becomes a significant factor influencing the overall reservoir performance [81,176].
As shown in Fig. 6 in single porosity models, hydrocarbon flow only occurs between matrix blocks neglecting fluid flow in reservoir fractures. However, in the dual-porosity approach, fluid flows to the well through fractures from the matrix blocks as they possess a much less permeability magnitude than the fractures. This indicates that the fluid flow can be exchanged between fractures and matrix blocks. Furthermore, in dual permeability models, fluid flow occurs through fractures, through the conductive matrix blocks, and between fracture and matrix cells. Accounting for fluid flow in these paths adds more complexity to fluid flow calculations and computational cost. The discrete fracture modeling approach accounts for fluid flow in fractures only; discarding overflow occurs through matrix blocks and isolated fractures. Fig. 6 illustrates the fluid flow mechanism through matrix and fracture systems for the models mentioned earlier.
The main difference between dual continuum modeling and discrete fracture modeling is in the early time flow behavior as well as the matrix fracture fluid flow exchange. The change in saturation is characterized by different accuracy and efficiency levels [130].

Implementation of Dynamic Reservoir Models in Fluid Flow Characterization Studies
The Warren and Root modeling approach has been subjected to several validation and analysis studies. For example, Lu et al. [19] have conducted a comparison study of the Warren and Root model with a mathematical model for dual-permeability (pressure transient model) to evaluate a vertical well's performance. Traditionally, Warren and Root model involved a critical assumption that the flow between matrix and fracture is in pseudo-steady state conditions; this indicates a uniform pressure decline through the matrix blocks. However, it is an oversimplified assumption when there is a contrast between fracture and matrix systems' permeability. Reference [177] validated pressure transient data using two new techniques known as simulation and inversion techniques. In order to model the fluid flow behavior and pressure performance of a fractured reservoir, the boundary element method (BEM) with random vertical fracture distribution is used. In addition, an analytical solution has been introduced by Kuchuk et al. [178] to understand the pressure performance of a continuous and discrete fractured reservoir. Pressure interpretation showed that Warren et al. [59] model is inadequate for analyzing a discrete fractured reservoir's pressure performance.
Lugumanov [45] has used a DPDP model to characterize fluid flow in fractured reservoirs. A fine grid hydrodynamic model was constructed to calculate the effective permeability. A multiplier was selected to Figure 6: Fluid movement through matrix and fracture systems for different fractured reservoirs modeling approaches simplify the process and obtain the effective permeability of the fractured medium. The multiplier was expressed as: where K ef is the effective permeability of the fractured medium, ψ is the clearance of the fractured medium in a given direction, È f is the porosity of the fractured medium, K f is the permeability of the fractured medium, and k is the variable parameter that makes the multiplier closer to the porosity (value closer to 1) or the clearance (value closer to 0). This study assumed that the size of the inhomogeneities was much larger than the molecular kinetic dimensions. Also, the size of the inhomogeneities was much smaller than the distances at which the macroscopic parameters change significantly. Under the standard conditions (at the surface), gas's solubility is assumed to be zero. Sensitivity analysis in numerical simulation indicated that the fractured reservoir's effective permeability is a function of the pore medium's absolute permeability.
Rueda Cordero et al. [111] stated that to appropriately characterize fluid flow in fractured reservoirs, integrating more than one numerical method became indispensable for planning new field development strategies and understanding the physical processes of fluid flow in NFRs. Besides, hybrid fluid flow modeling, such as integrating DFNM and dual continuum methods, is a reliable technique that provides accurate and realistic results. Lee et al. [179] have explicitly modeled long length fractures as 2D planes crossing multiple grid blocks. The calculation of transport parameters between fracture systems and discretized matrix blocks in this study was done systematically. Li et al. [129] have developed the technique proposed by Lee et al. [179] to circumvent all challenges involved with fracture object gridding. They introduced a different distinct model to simulate the flow of hydrocarbon in fractured reservoirs. To represent the matrix blocks, they used a conventional structured grid. They computed the intersection of fracture objects with the matrix blocks to introduce additional volumes of fracture control. The three main assumptions associated with these models are; perpendicular linear variation of pressure in the direction to each fracture, zero capillary pressure, and linear relative-permeability curves. DFNMs are found to be accurate as well as insensitive to grid orientation. Furthermore, the results were consistent and in good agreement with complex fracture patterns. In contrast, the conventional dual-permeability model indicated less satisfactory results; nonetheless, it showed a reasonable accuracy to model the flow behavior of a moderate fractured system. DPDP and DFNM are used in combination to represent fracture networks in finite element models [63]. Utilizing this hybrid technique may provide details about the reasons for early water breakthrough in fractured reservoirs, which, according to Bratton et al. [180], might be caused by high fracture density and permeability. Additionally, barriers can be formed by mineralized and cemented fractures, which prevent fracture networks' connectivity and result in poor fluid flow. This issue can only be detected using more than one numerical method [181]. Discrete fracture models are considered more realistic in representing a fractured system and view the explicit impact of individual fractures on hydrocarbon flow [13,104,105].
On the other hand, a more simplified way to characterize fluid flow in a fractured medium is to use DPDP continuum models, which are known by their distinct features like less running time and lower computational cost than DFNMs. In the case of simulating various fracture sizes in high-density fractured reservoirs, the DPDP models can help conduct the modeling process because the representation of each fracture is not needed [182]. Taabbodi et al. [81] stated that the dual-porosity technique was the most effective and popular approach to characterize hydrocarbon flow in the fractured medium before the recent rapid growth in computational power. In this study, fracture and matrix systems were counted as two separate continua with different properties. Furthermore, Bourbiaux et al. [160] stated that it is better to use hybrid modeling methods to characterize and evaluate fluid flow transport in NFRs. For example, combining continuum modeling with the DFNM approach [161], where large-scale fractures and major faults are modeled directly with DFNM, smaller fractures and vugs are described using a continuum modeling approach. These hybrid models are used to combine stochastic-based continuum modeling to quantify the uncertainty in specific geometrical properties of significant fractures and faults [18].
Many published studies in the literature, such as [173,[183][184][185][186], have indicated the role of geomechanics in fluid flow modeling and their impact on reservoir productivity. To capture fractured media's effect on hydrocarbon recovery in naturally fractured reservoirs, it is essential to couple both geomechanics and other fluid flow characterization techniques. The coupling process may come along with extra computational time and memory space. A more complex simulation model is needed since more factors are taken into consideration. The geo-mechanical part uses the equivalent continuum modeling approach that considers the deformed properties of the rock and fracture. Since most of the hydrocarbon is maintained within the matrix medium, the fracture porosity variation is deemed to be less critical for the coupled problem.
Luo et al. [187] proposed an integrated technique to describe the matrix and fracture influence on fluid flow. They coupled a parallel geomechanics model with a dual-porosity/dual permeability black oil model, as shown in Fig. 7. An iterative coupling process has been used to capture interactions between solid and flow. In order to handle large-scale fractures and more prolonged runtime problems, parallel computing is employed. Geo-mechanical effects on the reservoir pressure distribution were demonstrated using numerical experiments since different intensities of fractures have different directions and anisotropy behavior. In addition, to perform scalability behavior testing, it is required to use a large-scale model with very high grid cells number on multiple processors. Integrating more than one modeling approach can be useful and efficient in predicting and analyzing hydrocarbon production. The outcomes of the study also summed that this method could be used for full-field scale simulation studies. Furthermore, the results indicated an encouraging speedup of simulation time.
The pressure transient analysis can be associated with the DFNM approach for the permeability estimation of fractured reservoirs. Abdelazim et al. [71] have introduced a model to perform permeability calculations. Their methodology was as follows; first, the soft and hard data were analyzed to provide characteristic fracture details such as orientation, intensity, and fractures' dimension, after which fractures were distributed spatially using the geostatistical technique. Then fracture attributes were optimized through different realizations using gradient-based techniques. Finally, an innovative 3D discrete fracture  [187] model was built to estimate the pressure and its derivatives. A good match was obtained between the simulated pressure change and the pressure derivatives. It improved the modeling process by accurately capturing the complex subsurface fracture pattern. Operators can use the obtained details to design an efficient field development plan.
The approach of constructing a discrete fracture model has been widely illustrated in many studies, including [3,[114][115][116][117][118][119]. Based on different numerical simulation techniques, two discrete-fracture models have been presented by Moinfar et al. [13] to simulate the fluid flow behavior in naturally fractured reservoirs. The first model is represented by Unstructured Discrete-Fracture Model (USDFM), where grids around the fracture are characterized by unstructured cells with local refinement [127,128]. In contrast, fractures in the second model are incorporated in a conventional structured matrix grid using EDFM [188]. Both models were constructed to characterize a fractured medium's degree of heterogeneity in a more accurate way than traditional dual-permeability models. In both approaches, arbitrary fracture orientation was approximated by a horizontal rectangle plane that is orthogonal to the bedding. This approximation is justified because the fracture angle is very high concerning bedding [189]. This method might be used for fractures that possess arbitrary orientations in three dimensions.
In discrete models, a fractured medium can be separated into two parts; one describes the discrete fracture network, and another represents the permeable or impermeable matrix blocks. In addition, geometrical data such as fracture locations, shapes, and apertures are taken into consideration. Numerically, distinctive grid system designs are used to explicitly model fracture networks utilizing the fracture's width scale [190][191][192]. The grid generation of discrete fracture models is a crucial step that highly depends on the unstructured gridding technique. However, complex fracture network gridding may not be achieved using the available general-purpose grid processing tools. Therefore, Mallison et al. [193] have addressed this limitation by introducing a new algorithm for grid generation to model hydrocarbon flow in a fractured medium. Mallison et al. [193] raised an issue associated with selecting the precise position and geometry of fractures in NFRs, which is considered as an inherent uncertainty. Because of this issue, in their study, geological and flow simulation models were constructed without defining an exact fracture geometry. Grid generation algorithms and specified grid resolution are used to capture the geometric features, maintain a good quality of the grids, and honor fracture geometry. This approach helps to build multiple simulation models for a single characterization study.

Governing Equations for Fluid Flow in Naturally Fractured Reservoirs
The motion of fluid in NFRs is generally governed by the same fundamental principles of classical mechanics and thermodynamics based on momentum, mass conservation, and energy equations. Nonetheless, in conventional reservoirs, Darcy's law is employed instead of the momentum equation [56]. The following section will highlight the fundamental mathematical equations and their assumptions.

Continuity Equation
According to Barros-Galvis et al. [9], to model fluid flow in fractured reservoirs, mass conservation law (or continuity equation) can be used. This equation relates the mass rate change per unit volume in pore space with the sum of the fluid at the inlet and outlet points to provide the total mass accumulation in the control volume. The continuity equation can be expressed using a derivative or integral equation, which can be written as: where È f is the fracture porosity, r is the divergence operator, ρ is the fluid density, and ν is Darcy velocity.
Sarkar et al. [56] and Zimmerman et al. [194] argued that when the fluid is incompressible (the conservation of mass and volume are equivalent) and the fractured system lacks any source/sink, the continuity equation can be expressed as: where V is the velocity vector.
Eqs. (2) and (3) can provide a modeling description for the fluid movement through fractures. This set of equations might be used in simple flow cases, such as flow approximated by the parallel plate model for a single isolated fracture. However, in real reservoir cases, i.e., when fracture geometry deviates from the idealization, these equations might be solved using numerical solutions. Numerical models are constructed to simulate fluid flow through fracture networks and obtain solutions to pressure and velocity variations with space and time.

Equation of Motion
The general description of fluid dynamics in fractured reservoirs is governed by the Navier-Stokes equations based on fundamental physical laws of mass, momentum, and energy conservation at the microscopic level over the fracture void space. Under typical reservoir conditions, it is assumed that reservoir fluids are incompressible and Newtonian. In addition, it is assumed that the fluid density and viscosity are constant under steady-state conditions, and fracture walls are impermeable. The equation of motion is given by the Navier-Stokes equation, which can be written as [53,195,196]: where rP is the pressure of the moving fluid at each point, μ is the dynamic (or absolute) viscosity of the fluid, r 2 is the Laplacian operator and f e is the external force per unit volume acting on the flowing fluid.
Hydrocarbon flow in fractured reservoirs, particularly in narrow fractures, might be considered laminar, assuming that the fluid flows at a relatively low speed. Hence the effect of inertial forces could be neglected when solving the equation of motion. The effect of inertial and viscous forces can be determined by the dimensionless Reynolds number (Re), which can be stated as: where R e is Reynolds number and D À p is the average pore diameter. For natural fractures, the Reynolds number and the basic Darcy equation might be integrated to form Eq. (6) as: where q is the oil flow rate, a is the fracture aperture, A is the area, and Φ is the total porosity.
where k is the total permeability, and Φ is the flow potential.
A smaller Reynolds number indicates that the viscous effect resulting from the shear motion dominates the inertial forces associated with acceleration or deceleration of fluid particles in the fractured networks. This implies that the inertial terms, ρ(V⋅∇)V describing the fluid motion, can be omitted from Eq. (4).
The governing equation of motion will then be known as Stoke's equation and can describe what is called creeping flow. The time-dependent Stoke's equation might be written as: where: η = μ/q is the kinematic viscosity of the fluid.

Darcy's Law
According to Barros-Galvis et al. [9], Lu et al. [19], and Luo et al. [187], Darcy's law can be utilized to describe fluid flow through a fractured medium. It is commonly used to model the inner boundary condition for the corresponding radial flow model. In fractured reservoirs, Darcy's law can also be used to represent the relative fluid velocity using the following equation: where k i is the absolute permeability, k i ra is the relative phase permeability, m i a is the phase viscosity, B i a is the formation volume factor, p i a is the phase pressure, q i a is the fluid density, g is the gravitational constant, and H is the reservoir Depth.
Fluid flow in fractured reservoirs can also be represented using cubic law, which is given by: where w is the fracture aperture, p 1 is the outlet pressure, p 0 is the inlet pressure and L is the core length.
Based on Darcy's Law, the flow rate in fractured reservoirs can also be expressed as: where k x is the effective permeability, A is a constant.
Comparing Eq. (10) with Eq. (11), the intrinsic permeability for a single-phase flow in a single fracture plane is given by: Let H be the height of a grid block, and the effective permeability can be expressed as: 7.4 Partial Differential Equations for Two-Phase Flow (PDEs) De Souza et al. [197] argued that partial differential equations for two-phase fluid flow might be formed using a set of fundamental governing equations. These equations may include mass conservation, an expression for momentum balance, and equation of state. To obtain the pressure equation for two-phase flow, several assumptions must be considered, such as a heterogeneous reservoir, constant permeabilities, low or constant rock compressibility, dry gas reservoir, and Newtonian fluid. The mass conservation equation can be expressed as: where ν is the superficial fluid velocity, q m is the source term, v b is the bulk volume.
Tiab et al. [198] stated that viscosity must be modeled as a function of pressure and temperature for gas flow. They used Darcy's law for the low-velocity condition to express the conservation of momentum as follows: where k is the permeability, g is the gravity acceleration, and rH is the depth gradient.
The proposed dual-porosity model by Warren et al. [59] represents the fluid flow with two different mediums [199]. However, the primary concern consists of formulating the matrix-fracture transfer functions. Under the pseudo-steady-state condition, Warren et al. [59] introduced an equation for the overall matrix-fracture fracture flux as follows: where p m is the average matrix pressure, p f is the average fracture pressure, k m is the matrix permeability, and σ is the shape factor.
The dual-porosity reservoir comprises of two porosity systems, the primary matrix porosity and the secondary fracture porosity [84,85]. In order to analyze fluid flow in the matrix and fracture systems, it is crucial to use two main dynamic parameters, which are known as fracture conductivity, λ, and fracture storativity, ω. The fracture storativity measures the percentage of fracture contribution to reservoir fluid storage. A smaller value of ω indicates a more fracture effect on pressure variations. However, when ω value is close to unity, the reservoir behaves like homogeneous that made up solely of fractures. Fracture conductivity measures the ability of the matrix blocks to allow fluid to flow into the fracture networks. λ is related to the contrast in fracture and matrix permeabilities, i.e., the scale of heterogeneity. A smaller value of the fracture conductivity indicates a delay in the total system flow (transition period) [200,201]. Warren et al. [59] stated that two parameters are adequate to model a fractured reservoir's behavior (v is a dimensionless storage coefficient, whereas λ is the inner porosity flow parameter). These parameters can be expressed as: The fracture storativity ratio can be expressed as: where v is a dimensionless storage coefficient, Φ is the porosity, c t is the total compressibility, f denotes to fracture, and m denotes to matrix.
The Fracture conductivity is given by: where λ is the inter-porosity flow parameter, σ is the shape factor, r w is the wellbore radius, k e is the effective permeability of the fractured system.
The shape factor in Eq. (18) was subjected to many pressure transient studies and generally can be defined as an exchange coefficient to express matrix-fracture transfer. Warren et al. [59] have provided a geometrical definition of the shape factor as a term that reflects the matrix block's shape and represents the fluid flow between the two continua in fractured reservoirs. Their shape factor was expressed as: where n is the number of standard fracture sets, and l is the characteristic dimension of the heterogeneous region, which can be estimated from the surface-volume ratio: For n = 2 For n = 1 l ¼ 1 where a; b and c are the distance between fractures.
Kazemi et al. [2] defined the shape factor for three-dimensional NFR systems as follows: where L is the width of the matrix block in the three directions. Eq. (23) was based on the assumption that the pressure and saturation are uniforms within the girdblock. Coats [202] has defined the shape factor based on a pseudo-steady-state assumption as follows: Lim et al. [203] developed a shape factor for an isotropic rectangular system similar to Coast after reaching the pseudo-steady state condition. Their expression was defined as: The shape factor is considered as the most relevant feature in a dual-porosity model to account for matrix fracture flow. It is described mathematically by a transfer term. This term is highly dependent on matrix block shapes, flow regime (e.g., transient or pseudo-steady state), fracture pressure depletion, and production mechanism (e.g., gravity drainage, capillary imbibition, and capillary continuity). These factors are represented by one value known as the shape factor [204].

Conclusions
This paper has reviewed the modeling approaches which were introduced over the last decades to characterize fluid flow in naturally occurring fractured reservoirs; the following is concluded: Previous techniques present difficulties in modeling NFRs efficiently and accurately, particularly with multiple complex, large-scale fractures. Single continuum models are unable to effectively characterize the internal complexity of NFRs because of irregular reservoir shapes and different fracture scales. Therefore, a large number of theoretical methods and new technologies have been introduced to improve the modeling processes. Dual continuum models can better simulate fluid flow in a fractured porous medium. They are more common in the oil and gas industry due to their simplicity and lesser computational requirements than DFNMs. However, the fracture scale may influence the modeling process by adding more complexity. In comparison, discrete fracture modeling methods have been applied in some naturally fractured reservoir characterization problems. This approach seems to be an appropriate method for fluid flow characterization in reservoirs that contain a smaller number of fractures. However, it is challenging to run multiphase and full-field simulation studies due to many factors such as reservoir complexity, a large number of fracture sets, and flow paths, requiring extensive computational and memory power. Hybrid modeling is an exciting avenue for future research to produce a more accurate, less uncertain, and efficient reservoir model. The criterial selection of an appropriate reservoir model for fluid flow characterization in NFRs is highly dependent on the purpose of the simulation study, i.e., whether to describe fractures explicitly by a discrete fracture model or characterize them implicitly by an effective continuum model. However, it is not a binary choice, as there are several methods available to be applied solely or in combination.